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NUMERICAL  SOLUTION  OF  PLANE  VISCOUS  SHOCK  REFLECTIONS 

Tzu-Sien  Shao,  Ph.D. 
Department  of  Mathematics 
University  of  Illinois,  19^5 

ABSTRACT 

A  numerical  method  is  presented  for  the  solution  of  two-dimensional, 
viscous  fluid  flow  problems .   Emphasis  is  placed  on  the  formulation  and  solution 
of  weak  shock  reflections.   The  basic  equations  used  are  the  time -dependent 
Navier-Stokes  equations  without  heat  conduction.  A  version  of  the  shock  re- 
flection theory  dividing  flow  into  various  regions  is  discussed.  A  "behavior 
similar  to  experimental  results  is  displayed  "by  the  numerical  results  of  this 
theory,  although  they  do  not  quite  agree  with  the  experimental  data.  To  show 
that  the  Navier-Stokes  description  of  the  fluid  flow  is  an  adequate  model  for 
the  weak  shock  reflection  problems,  equations  of  motion  are  numerically  inte- 
grated by  a  difference  scheme  to  obtain  asymptotically  stationary  solution.  A 
transformation  of  independent  variables  is  first  made  to  reduce  computation 
time.   Stability  analysis  is  performed  to  give  an  upper  bound  for  the  time 
step.   Numerical  results  are  given  for  asymptotically  stationary  weak  shock 
reflections.   Reflected  shock  angles  for  both  regular  and  Mach  reflections 
obtained  agree  well  with  those  given  by  shock  tube  experiments. 


CHAPTER  I 
INTRODUCTION 

It  is  known  that  almost  all  hydrodynamical  phenomena  cause  a 
development  of  "discontinuities/1  so-called  shocks  and  boundary  layers,  sooner 
or  later.   The  essential  difference  between  a  shock  and  a  boundary  layer  may 
be  picturesquely  described  as  follows.  A  streamline  that  enters  a  shock  comes 
out  again  downstream,  whereas  a  streamline  that  enters  a  boundary  layer  remains 
within  a  region  of  shear.   Throughout  this  work,  we  shall  consider  shocks  only. 
Boundary  layers  will  be  ignored  completely.   Physically  speaking,  the  shocks 
are  regions  in  the  flow  where  higher  order  derivatives  of  certain  flow  varia- 
bles containing  small  coefficients,  e.g.,  coefficients  of  viscosity  and  heat 
conduction,  play  an  important  role.   Due  to  mathematical  difficulties,  these 
derivatives  are  usually  neglected  to  simplify  the  partial  differential  equations 
describing  the  flow.   This  causes  the  flow  variables  to  become  discontinuous 
across  the  shocks  which  now  are  curves  of  zero  thickness.   Using  the  integral 
form  of  the  equations  of  motion,  the  Rankine-Hugoniot  equations  can  be  derived 
to  determine  the  jump  conditions.   This  is  the  classical  approach.   It  works 
quite  well  for  simple  shock  problems. 

When  a  plane  shock  strikes  a  smooth  rigid  wall,  it  gives  rise  to 
regular  reflection  for  small  angles  of  incidence  and  to  an  irregular  type  of 
reflection,  known  as  the  Mach  reflection,  for  large  angles  of  incidence.   In 
regular  reflection  or  the  two-shock  case,  the  reflected  and  the  incident 
shock  meet  at  a  point  on  the  wall  (Fig.  l.l).   The  flow  direction  changes 
through  the  incident  and  reflected  shocks  must  be  equal  and  opposite  to  satisfy 
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the  conditions  of  flow  parallel  to  the  wall.   The  Mach  reflection  is  charac- 
terized by  the  intersection  of  three  shocks  at  a  point  detached  from  the  wall 
(Fig.  l.l).   This  is  known  as  the  three-shock  case.   The  shocks  are  called 
incident,  reflected  and  Mach  shocks  respectively. 


REGULAR  REFLECTION 


MACH  REFLECTION 


INCIDENT 
SHOCK 


^REFLECTED 
SHOCK 


////////J////////// 

WALL 
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REFLECTED 
SHOCK 
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FIGURE  1.1   REGULAR  AMD  MACH  REFLECTIONS 


For  this  class  of  shock  reflection  problems,  theories  based  on  the 
classical  approach  do  not  work  very  well.   In  particular,  they  fail  badly  for 
Mach  reflection  with  weak  incident  shocks.   There  is  a  wide  discrepancy 
between  the  reflected  shock  angles  predicted  by  the  three-shock  theory  and 
Smith's  [  1^+]  experimental  results  obtained  in  shock  tubes.  Another  anomaly 
'irst  observed  by  Smith  is  the  persistence  of  regular  reflection  somewhat 

nd  the  extreme  angle  where  the  two-shock  equations  cease  to  apply.   Since 
three-shock  theories  deal  with  shocks  in  the  neighborhood  of  the  inter- 
Lnt,  we  shall  refer  to  them  as  "local"  theories.   Presumably  local 
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theories  do  not  work  because  they  do  not  consider  the  influence  of  the  wall 
and  the  flow  downstream  of  the  intersection  region.   One  then  hopes  to  include 
these  effects  by  enlarging  the  intersection  region.   Theories  based  on  this 
idea  will  be  referred  to  as  "semi-global"  theories.   The  flow  field  can  be 
divided  into  various  domains  to  be  treated  separately  with  matching  conditions 
along  the  boundaries .  A  version  of  the  semi-global  theory  based  on  the  Navier- 
Stokes  description  of  fluid  flow  will  be  presented  in  Section  2.2.   It  does 
not  yet  give  satisfactory  results,  although  it  seems  to  work  better  than  local 
theories.   Therefore  one  must  solve  the  whole  problem  including  boundaries  and 
downstream  flow  flied.   This  gives  rise  to  a  "global"  approach. 

The  partial  differential  equations  of  compressible,  viscous  and 
heat  conductive  fluid  flow,  namely  the  Navier-Stokes  equations,  with  boundary 
conditions  are  too  complex  to  be  solved  analytically.  Numerical  solution 
seems  to  be  in  order.   Many  authors  [3>8]  have  numerically  solved  stationary 
fluid  flow  problems  with  artificial  viscosity  and  no  heat  conduction  from  a 
global  point  of  view.   Due  to  numerical  difficulties,  they  have  integrated 
time-dependent  equations  of  motion  to  obtain  asymptotically  stationary 
solution.   For  simple  shock  problems,  their  results  agree  quite  well  with 
classical  theory  and  experimental  data.   However,  no  method  seems  to  work  for 
three-shock  problems  with  weak  incident  shocks «   Since  the  true  viscosity 
effect  in  the  flow  is  more  important  in  this  case,  we  shall  include  it  in  the 
differential  equations  describing  the  flow.  A  numerical  method  for  solving 
asymptotically  stationary,  two-dimensional,  viscous  fluid  flow  problems 
from  global  point  of  view  will  be  presented  in  this  work.   Emphasis  will  be 
placed  on  the  formulation  and  solution  of  weak  shock  reflections.   The  purpose 
is  to  determine  whether  or  not  the  Navier-Stokes  equations  are  an  adequate 
model  in  this  case.   Since  the  shocks  are  weak,  the  corresponding  changes 
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in  flow  variables  are  small  throughout  the  flow  region.   Therefore  we  will 
neglect  heat  conduction  terms  in  the  equations  of  motion  without  effecting 
the  physics  of  weak  shock  reflection  problems  too  much.   The  coefficient 
of  viscosity  may  be  assumed  constant  for  the  same  reasons.   The  viscous  time- 
dependent  equations  are  nonlinear  and  parabolic,  and  therefore  very  time- 
consuming  to  solve  numerically  even  for  simple  shock  configurations.   Some 
efforts  have  been  made  to  make  our  method  more  efficient  and  easier  to  apply 
than  the  existing  ones.   The  main  features  of  the  present  methods  are:  [l)      it 
gives  a  faster  technique  for  numerical  solution  of  shock  reflection  problems 
and  it  is  directly  applicable  to  stationary  viscous  fluid  flow  problems  in 
general;  (2)  the  results  for  weak  shock  reflections  agree  very  well  with 
Smith's  shock  tube  experimental  data.   The  latter  is  a  strong  evidence  for  the 
adequacy  of  Navier-Stokes  description  of  shock  reflections., 

In  Chapter  II,  a  survey  of  previous  works  on  shock  reflection 
problems  will  be  presented  together  with  a  discussion  of  their  failure  for 
weak  three-shock  cases.  A  new  semi-global  theory  will  also  be  given  there. 
In  Chapter  III,  the  time-dependent,  non-heat-conductive  Navier-Stokes  equations 
will  be  transformed  into  a  streamline  coordinate  system.  A  more  detailed 
explanation  of  using  time-dependent  equations  to  obtain  asymptotically 
stationary  solutions  will  be  given.   The  equations  of  motion  will  be  written 
a  form  suitable  for  numerical  solution.   Chapter  IV  forms  the  main  text 
this  work.   Complete  numerical  procedures  for  solving  one-dimensional  and 
two-dimensional  viscous  fluid  flow  problems  will  be  derived  and  stability 
for  both  cases  will  be  performed.   Numerical  results  of  the  one- 
nal  viscous  shock  problem  and  plane  shock  reflections  made  on  the 
iversity  of  Illinois  ILLIAC  II  computer  will  be  reported.   Plots  of  flow 
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variables  drawn  "by  a  CalComp  67O  digital  plotter  will  also  be  shown.   Graphs 
consisting  of  data  obtained  from  local  theory,  semi-global  theory  and  our 
numerical  solutions  as  compared  to  experiments  will  be  given.   Good  agreement 
for  weak  shock  reflections  between  our  results  and  Smith's  experimental  data 
will  be  demonstrated. 


CHAPTER  II 
DEVELOPMENT  OF  SHOCK  REFLECTION  THEORIES 

2.1.   Failure  of  Local  Theories 

In  the  absence  of  a  complete  shock  reflection  solution  due  to 
mathematical  difficulties,  von  Neumann  [18]  suggested  a  local  three-shock 
theory  to  apply  to  the  neighborhood  of  the  intersection  of  the  incident , 
reflected  and  Mach  shocks.  The  idea  is  that  the  streamline  passing  through 
the  triple  point  divides  the  flow  into  two  contiguous  regions .   (This  dividing 
streamline  is  usually  referred  to  as  the  slip  stream.)   Combinations  of  Mach 
and  reflected  shocks  can  be  found  so  as  to  produce  the  same  static  pressure 
and  flow  direction  in  the  two  regions  downstream  of  the  intersection  point. 
Essentially  the  same  ideas  give  a  solution  for  regular  reflection  (two-shock 
theory)  which  corresponds  with  the  experimental  observations.   Unfortunately 
for  weak  incident  shocks,  there  is  an  apparent  wide  discrepancy  between  the 
reflected  shock  angles  predicted  by  the  three-shock  theory  and  Smith's  [Ik] 
experimental  results  obtained  in  shock  tubes. 

Stimulated  by  the  discrepancy  of  the  local  three-shock  theory,  a 
number  of  asymptotic  solutions  under  valid  limiting  conditions  have  been 
found  [1,9]-   However,  none  of  these  solutions  is  suitable  for  explaining  the 
failure  of  the  theory  for  reflected  shocks  with  finite  strength.  As  a  result 
of  the  limiting  conditions,  the  reflected  shock  is  forced  to  have  zero 
strength  at  the  triple  point. 

Later,  it  was  felt  [2]  that  perhaps  the  theory  fails  for  curved  shocks 
rather  than  weak  shocks.   Taub  [16]  worked  out  the  effects  of  curvature  and  its 

-6- 


-7- 

derivatives  in  the  neighborhood  of  the  triple  point.   Clutterham  [k]    calculated 
reflected  shock  angles  as  a  function  of  the  angle  of  incidence  and  the  incident 
shock  strength  based  on  the  curved  shock  theory.   Nevertheless  this  local 
theory  of  curved  shocks  still  suffers  from  any  weakness  that  the  three-shock 
theory  suffers.   This  is  due  to  the  fact  that  it  contains  the  three-shock 
theory  in  the  limit  when  the  neighborhood  under  consideration  shrinks  to  the 
triple  point. 

In  1959*  Sternberg  [15]  examined  the  effect  of  shock  thickness  on 
the  boundary  conditions  at  the  triple  point.   He  showed  that  at  the  intersection 
there  must  be  a  non-Rankine-Hugoniot  shock  zone  separating  the  three  R-H  shocks. 
According  to  him,  the  three-shock  theory  of  von  Neumann  and  Guderley  [7]  fails 
because  they  used  nonviscous  models.   With  viscosity  included,  there  is  a 
smooth  transition  of  flow  variables  inside  the  non-R-H  zone.   This  non-R-H 
zone  acts  as  a  buffer  zone  separating  R-H  shocks  in  the  upper  and  lower  domains. 
It  then  completely  eliminates  the  singularity  at  the  triple  point  appearing  in 
nonviscous  models.   However,  the  two-dimensional  character  of  the  non-R-H  zone 
introduces  a  theoretical  complication.   The  unknown  downstream  elliptic  field 
has  a  significant  influence  on  the  structure  of  the  shock  zone.   This  makes 
the  mathematical  problem  indeterminate  from  a  local  point  of  view. 

Later  Gear  [6]  tried  to  introduce  viscosity  into  the  curved  shock 

theory.   For  very  weak  shocks,  the  Mach  number  M  ahead  of  the  incident  shock 

1  2 
is  close  to  unity.   Writing  M  =  1  +  —  e  ,  he  expanded  flow  variables  in 

power  series  of  the  small  parameter  e.   For  the  structure  of  shock  regions, 

he  used  a  linear  approximation.   Retaining  only  low  order  terms  in  e   under 

consistent  assumptions,  a  system  of  nonlinear  algebraic  equations  was  derived 

relating  the  reflected  shock  angles  and  shock  curvatures.   The  only  numerical 
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solution  of  the  algebraic  system  found  prescribed  a  negative  shock  in  violation 
of  the  second  law  of  thermodynamics.   Sukurai  [13]  also  tried  to  use  power 
series  expansion  in  terms  of  radial  component  of  the  independent  variable  to 
integrate  the  viscous  differential  equations.   His  choice  of  boundary  conditions, 
however,  remains  to  be  justified. 

2.2,  A  Semi-Global  Shock  Reflection  Theory 

Since  all  local  theories  fail  for  weak  Mach  reflection  problems,  we 
propose  in  this  section  a  semi-global  theory  to  improve  the  results.   This 
theory  includes  the  effects  of  the  wall  and  the  slip  stream  behind  the  inter- 
section region.   The  shocks  are  treated  as  domains  with  finite  thickness  to 
avoid  singularities.  As  in  Gear's  work  [6],  the  flow  field  is  divided  into 
several  regions.   The  intersection  region  is  thought  of  as  a  mixture  of  incident, 
reflected  and  Mach  shocks,  outside  of  which  viscosity  can  be  ignored.   The 
boundary  condition  at  the  wall  to  which  the  Mach  shock  is  attached  will  be 
assumed  to  correspond  to  the  condition  that  the  wall  exerts  no  viscous  drag 
on  the  fluid.   This  condition  can  be  realized  when  two  identical  shocks  inter- 
sect each  other  symmetrically.   The  central  dividing  streamline  is  physically 
equivalent  to  the  wall  exerting  no  viscous  drag.   The  idea  is  to  integrate 
equations  of  motion  along  the  boundaries  of  the  intersection  region.  A  system 
of  nonlinear  algebraic  equations  relating  the  flow  variables  and  geometric 
variables  will  then  be  obtained.  Appropriate  solution  of  these  equations 
will  yield  reflected  shock  angle  and  Mach  shock  height  for  given  incident 
shock.   For  simplicity,  we  use  the  following  two-dimensional,  stationary  and 
non-viscous  equations  of  motion  for  a  perfect  fluid  in  rectangular 
coordinates  [17]. 


(l)   Continuity  equation 


c*  (pu)  +  "oy  ^  =  °  ; 


(2.2.1) 


(2)     Momentum  equations 


(pu   ,+  p>  +>x-   (puv)    =  0    , 


bx 


% 


(2.2.2) 


3x"  (pUv)    +  ¥  (pV     +  P)    =        j 


(2.2.3) 


(3)      Energy  equation 


pu( 


2      2 
u  +v 


+  -2-£) 

7-1   p-J 


fr 


2      2^' 
/U  +v         _7     V\ 

pv( — - —  +  — —  — ) 
w    K      2  7-1  p\ 


=  0  (2.2. k) 


where  7  is  the  ratio  of  specific  heats .  We  shall  discuss  the  semi-global 
theory  in  two  cases  using  different  geometric  models  for  the  intersection 
region.   They  are  as  follows. 

(A)   1-Region  Model 


■8X— H«s—  g 
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rrrrrm 
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FIGURE   2.1        SEMI-GLOBAL  1-REGION  MODEL 
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Here  we  consider  a  model  (Fig.  2.l)  where  the  intersection  region 
is  connected  to  the  incident  and  the  reflected  shocks  and  extended  to  the 
wall.   The  velocity  component  normal  to  the  wall  at  A'  and  D  is  taken  to  "be 
zero.   The  shocks  are  assumed  to  be  stationary  and  oblique  to  the  flow 
direction.   The  flow  conditions  in  zone  1  are  assumed  uniform  and  given.   By 
specifying  the  incident  shock  strength   £T  and  shock  angle  CC  we  can  calcu- 
late the  flow  conditions  in  zone  2  by  using  the  conditions  for  a  stationary- 
oblique  shock  [10].   Similarly,  flow  conditions  at  B  can  be  determined  if  the 

reflected  shock  angle  a     is  assumed  in  computation.   Thus  the  unknowns  in 

K 

this  model  are  the  flow  variables  at  D,  the  reflected  shock  angle  a     and  the 

K 

height  of  the  intersection  region  h.   There  are  actually  5  unknowns,  namely 

p  ,  un,  p_,  &     and  h  since  v  is  taken  to  be  zero.   If  we  integrate 
D        D        D        t\  D 

Eqs.  (2.2.1)  -  (2.2.  k)    over  the  region  ABDA.'  with  respect  to  x  and  y,  we  can 

get  k   equations  relating  the  above  unknowns.   We  shall  illustrate  the  inte- 
gration procedure  by  using  Eq.  (2.2.1)  as  follows. 

Integrating  Eq.  (2.2.1)  over  ABDA.'  with  respect  to  x  and  y,  we  have 


B 
(pu)dy 


B    fB    B 

+  /   [pv]   dx  -  0  (2.2.5) 

A   JA  U 


In  order  to  numerically  handle  the  above  integrals,  we  first  observe  that  the 

intersection  region  can  be  thought  of  as  having  a  mirror  image  with  respect 

to  the  wall.   (This  is  equivalent  to  assuming  that  the  wall  exerts  no  viscous 

drag.)  The  following  modified  Simpson's  rule  can  be  derived  to  take  advantage 

of  the  symmetry 
x 
f         f(x)dx  =  I  [f(xQ)  +  2f(x1)]  -  |-  f(lv)U)  ,  (2.2.6) 

o 
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where     h  =  x 

and       x  <  i   <  x 

o 


o    o 


1 


This  formula  can  be  used  for  integrations  from  D  to  B.   Integrals  from  A  to 
B  are  handled  by  the  standard  trapezoidal  rule.   Using  the  conditions 


v.  =  v . ,  =  v  ~  0,   Eq.  (2.2.5)  now  becomes 


§  [(pu)B  +  2(pu)D]^  +  -±   (pv)I  +  -f  [(pv)1   +  (pv)B]  =  0  ,     (2.2.7) 

where  5  and  6  are  the  oblique  thicknesses  of  the  incident  and  reflected 
I      K 

shocks  respectively.   They  can  be  calculated  by  using  the  solution  of  one- 
dimensional  viscous  shock  problem  [17]-   Further  simplification  of  Eq.  (2.2.7) 
using  the  uniformity  condition  of  the  incoming  flow  in  zone  1  leads  to 


—  (pu)D  =  h(pu)A  -  ^(pu)B  -  —  (pv)B 


6  +5 

I  R  /   x 
-g-  (pvJj.  • 


(2.2.8) 


Similarly,  Eqs .  (2.2.2)  -  (2.2. k)    can  be  integrated  and  written  in  the  form 


2h  /   2   s    _  ,   2   v    h  t      2   s    &R  t        v     I+5R  /    % 
—  (pu  +p)D  =  h(pu  +p)A  -  -  (pu  +p)B  -  —  (puv)B  -  — - —  (puv). 


I 
(2.2.9) 


5I+5R      h  (        ,  5I+5R  ,  2      s  5R  ,   2   ,    &R 

-2~  PD  -  3  (puv)B  +  -5—  (pv  +p)I  +  -  (pv  +P)B  -  -  pA 


(2.2.10) 
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respectively 


(2.2.11) 
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We  now  have  k   equations  in  5  unknowns.   The  extra  equation  can  be 
supplied  by  assimung  that  the  flow  behind  the  intersection  region  is  shear- 
less.   That  is  to  say,  the  horizontal  velocity  components  at  B  and  D  have 
the  same  value.   We  thus  have 

uD  =  uB  .  (2.2.12) 

Eqs.  (2.2.7)  _  (2.2.12)  form  a  system  of  5  nonlinear  algebraic  equations  in 

5  unknowns  p_,  u^  p  ,  GC     and  h  where  Ot     appears  implicitly  through  the  flow 
D        D      D        K  K 

conditions  at  B.   Further  manipulations  of  the  equations  can  reduce  them  to 

one  single  equation  involving  a  .  A  numerical  solution  was  found  using  the 

n 

binary  chopping  technique,  giving  a     for  given  %     and  a    .   The  Mach  shock 

height  h  can  also  be  calculated.,   The  results  for  £   =  0.9;  0.8,  0.7  and 

a     =  k6   (l  )70  are  given  in  Table  2.1.   The  shock  angles  are  also  plotted  in 

Figures  k.^   -  ^■•9-      Since  the  equations  are  very  nonlinear,  several  extraneous 

solutions  were  discarded.   Fortunately  the  extraneous  solutions  can  easily 

be  distinguished.  They  all  fall  into  the  following  categories:   (l)  a 

corresponds  to  no  reflected  shock,  i.e.,  uT  =  u_;  (2)  a   corresponds  to 

1     B     '   K 

negative  value  of  h;  and  (3)  0!  is  less  than  a    . 

n  I 

(B)   2-Region  Model. 

Encouraged  by  the  1-Region  model,  we  consider  a  more  complicated 
and  presumably  more  realistic  model.   There  is  a  slip  stream  extending  from 
the  intersection  region  through  BC .   (See  Fig.  2.2).  All  assumptions  and 
conditions  are  taken  to  be  identical  to  case  (A).   The  unknowns  here  are 
low  variables  at  C  and  D,  the  reflected  shock  angle  QL,,  the  height  of  slip 
tream  s  and  the  distance  from  slip  stream  to  the  wall  h.  With  v  taken  to 
be  zero,  we  have  10  unknowns,  namely,  pQ,   u  ,  v  ,  p  ,  pD,  u  ,  pD,  aR,  s  and  h. 
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FIGURE  2.2   SEMI-GLOBAL  2-REGION  MODEL 


As  in  case  (A),  we  integrate  Eqs .  (2»2.l)  -  (2.2.^)  over  regions  ABDA"  and 


A'CDA"  separately  in  x  and  y  directions.   Using  the  conditions  v.  =  v.  ,  = 


V  =  VD  =  °'   Eq°  C2-2'1)  Sives 


h 


2h 


.  (pu)c  +  —  (pu)D  +  -  (pu)c  =  h(pu)A  +  s(pu)A  -  -  (pu)B  -  —  (pv)I 


2"  [(PV).B  +  (p*)^ 


(2.2.13) 


and 


—  (pu)D  =  h(pu)A  -  -  (pu)c 


5I+5Rr   , 
-g—  (pv)c 


(2.2.1^) 


respectively  for  integrations  over  regions  ABDA"  and  A'CDA".   Substituting 
Eq.  (2.2.11+)  into  Eq.  (2.2.13)  and  simplifying,  we  can  eliminate  p  ,  u  ,  and 
h  and  set 
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5I+SR 


5I+6R 


|  (pu)c  -  -—^   (pv)c  =  s(pU)A  -  -  (pu)B —  (pv)I  -  —  (pv)B  . 

(2.2.15) 
Similarly,  Eqs.  (2.2.2)  -  (2.2.4)  can  "be  integrated  and  simplified  to  yield 
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(2.2.16) 
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(2.2.17) 
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(2.2.18) 


respectively.   In  Eqs .  (2.2.15)  -  (2.2.18),  unknown  flow  variables  p  ,  u  ,  p 
and  the  distance  h  have  "been  eliminated.   We  now  have  k   equations  in  6 


unknowns  p  ,  u  ,  v  ,  p  ,  a,     and  s  where  a     again  appears  implicitly  through 
L    U    U    0    K  R 

flow  conditions  at  B.   To  supply  two  more  equations  for  the  problem,  we 
assume  that  flow  variables  across  the  slip  stream  satisfying  the  classical 
conditions.   These  conditions  can  be  derived  by   integrating  Eqs.  (2.2.1)  - 
(2.2.4)  across  BC  and  letting  s  -*  0.   They  are  simply 


PC   =  PB 


111       Ib 
uc 


(2.2.19) 
(2.2.20) 
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Eqs .  (2.2.15)  -  (2.2.20)  now  form  a  system  of  6  equations  in  6  unknowns. 

Further  manipulations  and  simplifications  can  again  reduce  them  to  a  single 

equation  in  a    .   Numerical  solutions  are  given  in  Table  2.1  for  |   =  0.9, 

0.8,  0.7  and  a     =  h6°(lO)70°.      In  Figures  4.  7  -  k>9,   we  shall  plot  the  results 

of  2-region  model  as  compared  to  1-region  model,  Smith's  experiments,  global 

theory  and  local  theory.   The  extraneous  solutions  discarded  in  this  case 

fall  into  the  categories:  (l)  CC  corresponds  to  no  reflected  shock;  (2)  QL, 

R  R 

corresponds  to  negative:  values  of  s  or  h;  (3)  oc     is  less  than  a,     and  (h)  a 

K  In 

corresponds  to  unreasonably  large  values  of  s  and  h. 

The  calculated  shock  reflection  angles  in  both  cases  do  not  quite 
agree  with  Smith's  experimental  data  as  we  can  see  in  Figures  ^.7  -  ^■•9* 
They  seem  to  give  better  results  than  local  theories  for  weak  three -shock 
intersections,  however.   Strangely,  examination  of  Table  2.1  and  Figures  ^.7  - 
K.y,    shows  that  the  1-region  model  consistently  gives  better  results  than  the 
2-region  model,  although  the  heights  of  intersection  region  in  both  cases 
agree  fairly  well  as  shown  in  Table  2.1.   We  do  not  yet  have  a  good  explana- 
tion for  this  phenomenon.   The  assumption  of  shearless  flow  behind  the  inter- 
section region  used  in  1-region  model  does  not  seem  to  have  much  physical 
justification.   This  more  or  less  arbitrary  assumption  might  accidentally 
give  better  answers.   In  the  mean  time,  conditions (2.2.19)  -  (2.2.20)  used 
in  2-region  model  may  not  be  proper  when  viscosity  is  present.  A  better 
understanding  of  the  viscous  effect  along  the  boundaries  of  the  intersection 
region  is  needed  to  improve  the  results. 

It  is  conceivable  that  some  approach  similar  to  the  semi-global 
theory  discussed  in  this  section  could  be  useful  for  shock  intersection 
problems.   The  flow  field  could  be  divided  into  various  regions  where  certain 
flow  parameters  may  be  more  important  than  others  in  a  given  region. 
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Integrating  equations  of  motions  in  these  regions  separately  with  matching 
conditions  along  the  boundaries,  we  could  obtain  flow  conditions  behind  the 
intersection  region  and  geometric  variables  associated  with  the  regions. 
With  the  lack  of  success  of  local  theories  and  the  present  incomplete  semi- 
global  theory ,  we  proceed  to  solve  the  shock  intersection  problem  numerically 
from  a  global  point  of  view.   This  forms  the  main  text  of  Chapters  III  and  IV. 

2.3'   Global  Theory  with  True  Viscosity 

As  we  have  seen,  all  local  three-shock  theories  as  well  as  the 
proposed  semi-global  theory  seem  to  fail  in  one  way  or  the  other.   The  question 
naturally  arises  whether  or  not  there  is  a  global  theory  with  true  viscosity 
which  would  work  for  all  incident  shocks.   Since  the  effects  of  wall  and  the 
flow  behind  the  intersection  are  believed  to  be  important,  any  global  theory- 
must  consider  them  in  the  boundary  conditions  of  the  whole  problem.  Analytical 
solution  of  the  viscous  equations  of  motion  with  boundary  conditions  appears 
out  of  the  question.  A  complete  numerical  solution  by  means  of  high  speed 
digital  computers  seems  to  be  in  order. 

Numerically,  the  appearance  of  discontinuities  makes  it  difficult 
to  integrate  the  basic  equations  throughout  the  entire  flow  region  under 
consideration.   Von  Neumann  and  Richtmeyer  [19]  proposed  a  method  for  numerical 
calculation  of  hydrodynamic  shocks.   They  introduced  the  so-called  "artificial 
viscosity"  to  smear  out  the  shock  regions.   Thus  the  shock  is  not  considered 
as  an  interior  boundary,  but  as  a  part  of  the  solution,,   The  idea  is  to  replace 
the  true  viscous  terms  in  the  differential  equations  by  simpler  dissipative 
expressions.   Later  Harlow  [8]  proposed  the  particle-in-cell  (PIC)  method  for 
numerical  solution  of  fluid  dynamics  problems.   Rich  [11]  recently  used 
explicit  Landshoff  type  viscous  pressure  as  well  as  implicit  dissipation  terms 
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arising  from  the  "repartitioning"  of  the  energy  and  momentum  in  cells  for 
numerical  calculations  of  detached  shocks  with  success.   All  these  techniques 
were  mainly  used  for  relatively  simple  problems  nevertheless.   For  shock 
reflection  problems,  they  are  still  too  time-consuming  with  the  existing 
computing  equipment. 

In  196^-,  Bur  stein  [3]  numerically  integrated  the  time -dependent 
inviscid  equations  of  motion  by  the  Lax-Wendroff  difference  method.   He 
obtained  results  for  asymptotically  stationary  regular  and  Mach  reflections 
in  air.   The  Mach  reflection  calculations  for  very  strong  shocks  agree  with 
experimental  photographic  data  obtained  from  wind-tunnel  tests.   This  is  to 
be  expected  as  the  local  three-shock  theory  has  already  proven  satisfactory 
for  strong  shocks.   However,  the  results  for  Mach  reflection  of  weaker 
incident  shocks  were  quite  disappointing  as  compared  to  Smith's  experiments. 

It  seems,  after  reviewing  previous  attempts  at  the  problem,  that 
the  true  viscous  terms  play  a  vital  role  for  the  solution  of  weak  shock 
reflections.   Numerical  solution  of  the  viscous  equations  can  give  a  global 
view  of  the  whole  problem.   Due  to  the  slowness  of  the  asymptotic  development 
of  the  shock  intersection  region  as  well  as  the  complexity  of  the  true 
viscous  terms  in  the  equations,  a  more  efficient  numerical  method  than  the 
existing  ones  is  needed.   We  will,  in  the  following  chapters,  present  such  a 
global  approach  with  its  numerical  solution  to  show  that  the  Navier-Stokes 
equations  are  an  adequate  description  of  the  flow  under  these  conditions. 


CHAPTER  III 
DIFFERENTIAL  EQUATIONS  SUITABLE  FOR  NUMERICAL  SOLUTION 

5»1«   Time -Dependent  Equations  Versus  Stationary  Equations 

The  class  of  plane  shock  intersection  phenomena  under  consideration 
is  stationary.   That  is  to  say,  the  value  of  a  flow  variable  at  a  given  point 
in  a  fixed  coordinate  system  does  not  change  in  time.   It  is  then  natural  to 
use  stationary  viscous  equations  of  motion  for  numerical  integration. 
However,  in  applying  global  theory  to  shock  reflection  problems,  there  are 
numerical  difficulties.  A  different  approach  is  needed.   In  this  work,  we 
shall  integrate  time -dependent  viscous  equations  in  Lagrangian  form  with 
initial  and  boundary  conditions  to  obtain  asymptotically  stationary  solutions. 
It  is  known  that  stationary  fluid  flow  equations  are  nonlinear  even 
without  true  viscosity  terms.   There  are  no  convenient  numerical  methods  to 
handle  the  nonlinearities  that  can  be  guaranteed  to  converge.   The  situation 
becomes  worse  when  viscosity  is  introduced.   If  one  considers  the  time- 
dependent  equations,  the  time  derivatives  of  the  flow  variables  appear 
linearly.   One  can  therefore  difference  the  equations  to  iterate  in  time. 
The  iterative  process  can  be  continued  until  values  of  all  flow  variables 
between  consecutive  time  cycles  differ  by  a  negligible  amount.   Since  the 
time  derivatives  are  then  small,  the  flow  is  the  solution  of  a  set  of  equa- 
tions which  are  disturbed  from  the  original  stationary  equations  by  a  small 
amount.   In  the  limit,  the  solution  thus  obtained  is  stationary.   On  the 
other  hand,  if  the  problem  is  initially  ill-posed,  one  may  not  get  convergence 
This  fact  makes  it  plausible  to  use  time -dependent  equations  of  motion. 
There  are  in  general  two  basic  ways  of  writing  time-dependent 
equations;  they  are  called  Lagrangian  and  Eulerian.   In  the  Lagrangian  form, 
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the  coordinate  system  is  moving  with  the  fluid.   The  main  advantage  of 
Lagrangian  techniques  has  to  do  with  the  characteristics  and  the  Courant 
stability  condition.   The  equations  of  motion  are  also  simpler  in  form. 
The  accuracy  of  Lagrangian  methods  is,  however,  doubtful  when  applied  to 
a  system  after  distortions  of  the  fluid  have  occurred.   In  the  Eulerian 
view  point,  the  coordinate  system  is  fixed  relative  to  the  observer  and  the 
fluid  moves  through  the  mesh  of  coordinates.   Eulerian  methods  have  the 
advantage  of  applicability  to  problems  with  arbitrary  distortions  or  slippage 
of  the  fluid.   They  also  suffer  from  several  disadvantages  such  as  the  false 
diffusion  arising  from  forcing  the  meshs  to  be  homogeneous.   It  is  then 
desirable  to  combine  the  features  of  the  Lagrangian  and  Eulerian  methods  and 
thereby  eliminate  some  of  the  disadvantages.   The  technique  we  shall  adopt 
is  to  first  write  the  time -dependent  equations  of  motion  in  Lagrangian  form. 
An  iterative  numerical  procedure  is  then  derived  from  these  equations .   The 
corresponding  difference  equations  are  iterated  in  time  to  obtain 
asymptotically  stationary  solution.   However,  we  are  interested  in  the 
stationary  solution  from  Eulerian  point  of  view.  Any  solution  of  the 
Lagrangian  equations  of  motion  alone  will  not  be  sufficient.   This  problem 
can  be  solved  by  "repartitioning"  fluid  cells  or  "interpolating"  fluid  mesh 
points  depending  on  the  basic  numerical  subdivision  used  in  the  process. 
The  basic  idea  is  to  re-evaluate  the  flow  variables  in  the  fixed  Eulerian 
mesh  after  each  Lagrangian  time  cycle.   The  transport  effect  of  flow  variables 
can  thus  be  obtained.  A  more  detailed  discussion  of  the  "repartitioning" 
and  "interpolating"  processes  will  be  given  in  the  next  section. 

5.2.  The  Cell  Method  and  the  Point  Method 

In  most  numerical  techniques  for  solving  fluid  dynamics  problems, 
the  fluid  region  is  subdivided  into  small  cells  or  mesh  points .   Finite 
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difference  approximations  to  partial  differential  equations  discribing  the 
flow  are  derived  and  applied  to  each  cell  or  mesh  point.   They  will  be  referred 
to  as  the  cell  method  and  the  point  method  respectively. 

In  1957;  Harlow  and  his  colleagues  [8]  developed  the  particle-in-cell 
method,  abbreviated  PIC  method,  for  the  numerical  solution  of  problems  involving 
the  dynamics  of  compressible  fluids.   It  combines  the  features  of  the  Lagrangian 
and  the  Eulerian  methods  and  thereby  eliminates  some  of  the  disadvantages. 
There  are  two  computing  meshes;  one  is  Eulerian,  the  other  Lagrangian.   The 
domain  through  which  the  fluid  is  to  move  is  divided  into  a  finite  number  of 
cells o   This  is  the  Eulerian  mesh.   In  addition,  the  fluid  itself  is  represented 
by  a  finite  number  of  particles  which  move  through  the  Eulerian  mesh,  repre- 
senting the  motion  of  fluid.   This  is  the  Lagrangian  mesh.  Associated  with 
the  mesh  points  of  each  system  are  certain  variables  whose  history  the  calcu- 
lation develops.   Thus  for  each  Eulerian  cell  there  is  kept  the  velocity,  the 
internal  energy,  and  the  total  mass  of  each  kind  of  material.   For  the 
Lagrangian  mesh  of  particles,  individual  masses  and  positions  are  kept.   The 
differential  equations  of  motion,  with  transport  terms  neglected,  are  written 
in  finite  difference  form  relative  to  the  system  of  cells.   The  transport 
effect  is  obtained  later  by  allowing  the  particles  to  move  through  the 
Eulerian  mesh  according  to  the  velocities  in  the  neighboring  cells. 

Later,  Rich  [11]  modified  the  PIC  method  for  the  case  of  a  continuous 
fluid.   He  applied  the  integral  form  of  the  conservation  laws  to  the  whole 
volume  of  each  cell.  As  in  the  PIC  method,  the  calculation  of  the  configura- 
tion of  the  fluid  at  time  t   ,  from  that  at  time  t  may  be  thought  of  as 

n+1  n 

proceeding  in  three  phases.   In  the  first,  the  fluid  element  in  each  cell  is 
assumed  fixed  and  no  flow  across  the  cell  boundaries  is  allowed .   By  applying 
conservation  equations,  intermediate  values  of  the  velocity  and  specific 
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energy  in  each  cell  are  obtained.   The  second  phase  consists  of  the  calculation 
of  the  mass  flows  across  the  cell  boundaries.  Mass  crossing  a  boundary  is 
assumed  to  carry  with  it  the  intermediate  velocity  and  specific  energy  of  the 
cell  from  which  it  came.   In  the  third  phase ,    the  new  cell  densities  are 
computed;  and  by  application  of  conservation  of  energy  and  momentum,  the  new 
velocities  and  specific  energies  are  subsequently  obtained.   Rich's  results 
on  the  supersonic  flow  of  air  past  a  variety  of  cylindrically  symmetric  objects 
agree  quite  well  with  experiments. 

In  1964,  Burstein  [3]  differenced  the  inviscid  conservation  equations 
at  each  Eulerian  mesl  joint  by  the  Lax-Wendroff  method.   It  is  based  on  the 

Taylor  expansion  of  the  vector  function  w(x,y,t  -1-  At)  so  as  to  include  the 

1   2 
second  order  term  —  At  w,   .   Since  the  difference  scheme  is  made  conditionally 

stable,  with  the  addition  of  this  term,  a  linear  analysis  allows  the  computa- 
tion of  the  linear  stability  limit.   The  time  increment  At  and  space  increments 

2    2  l/2 
Ax  and  Ay  are  related  to  the  flow  Mach  number  M  =  (u  +  v  )  '  /c  (c  being  the 


local  sound  speed)  by 

At  =  At  <   (M^l)-1  (     j 

Ax     Ay   -  gl/2  c  ^•J-LJ 

This  numerical  scheme  was  tested  for  various  values  of  incident  shock  strengths 
For  regular  reflection,  the  results  agree  quite  well  with  exact  solutions. 
The  Mach  reflection  calculations  for  weak  incident  shocks  did  not  compare  well 
with  Smith's  experiments. 

Now,  we  have  seen  that  the  PIC  method  treats  each  cell  as  composed 
T  many  individual  particles,  each  of  them  is  free  to  move  from  one  cell  to 
the  other.   Rich  considers  each  cell  as  a  unit  with  continuously  distributed 

it.  These  will  be  referred  to  as  the  "cell  method."  The  equations  of 
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motion  are  integrated  with  respect  to  the  cell  first,  and  then  the  integrals 
are  approximated  by  suitable  quadrature  formulas .   As  in  Burstein's  work,  we 
shall  use  the  "point  method"  for  numerical  integration.   The  flow  region  will 
be  divided  by  a  mesh.   The  Lagragian  equations  of  motion  will  be  differenced 
at  each  mesh  point  so  that  one  time  step  serves  to  displace  the  Lagrangian 
mesh  points  relative  to  the  Eulerian  mesh  points.   The  physical  positions  of 
the  displaced  mesh  points  relative  to  the  Eulerian  mesh  points  can  be  calcu- 
lated from  the  velocity  components «   Then  the  new  flow  variables  on  the  fixed 
Eulerian  mesh  can  be  evaluated  from  the  intermediate  values  of  the  flow 
variables  computed  in  Lagrangian  time  cycle  by  interpolation. 

3«3°   Differential  Equations  in  Lagrangian  Form 

The  theory  of  fluid  flow  that  we  will  use  is  based  on  the  Newtonian 
mechanics  of  a  small  body.   For  continuously  distributed  fluid  particles  in 
the  flow,  the  essential  part  of  Newton's  principles  can  be  described  by  the 
conservation  laws  of  mass,  momentum  and  energy  of  a  fluid  element  [17] •   For 
simplicity,  we  assume  the  fluid  is  not  heat  conductive.   From  the  result  of 
one-dimensional  flow  [16],  principle  features  of  the  problem  can  be  found 
under  the  assumption  of  no  heat  conduction.   For  the  viscosity  terms,  we 
assume  that  the  stress  tensor  is  a  linear  function  of  the  rate  of  strain 
tensor.   The  coefficient  of  viscosity  is  taken  to  be  constant  and  the  bulk 
viscosity  coefficient  is  assumed  to  be  zero=   These  assumptions  are  particu- 
larly justifiable  when  weak  shock  intersection  problems  are  considered.   Under 
these  assumptions,  we  can  write  the  time-dependent,  two-dimensional,  viscous 
equations  in  Lagrangian  form  [10]  as  follows. 

(l)  Continuity  Equation 

p  +  pu.  .  =  0  ,        i  =  1,2  (3-2.1) 

1,1 
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(2)  Momentum  Equations 


pu.  = 


u(u .  .+u  .  .  )  -  —  Liu.  .   o  .  .  -p6  .  . 


,   j  =  1,2     (30-2) 


where  u.  is  the  coefficient  of  viscosity  and  5.  .  is  the  Kronecker  delta, 
(3)   Energy  Equation 


p(e  +  g  u1ui) 


LLU.(u.   .+U.   .  )  -  —  liU,   ,U.-pu. 

J  i,J  J,i    3   k,k  1  ^  ij 


■p(e  +  -  u.u.  )u.  . 


(3.3.3) 


where  e  is  the  specific  internal  energy 


In  the  above,  we  have  used  a  dot  over  a  function  to  denote,  not  partial  dif- 
ferentiation with  respect  to  time  t  at  constant  (x  ,x  ),  but  rather 
differentiation  for  a  given  particle  whose  positions  changes  according  to  its 
own  velocity  vector.   This  is  known  as  the  Euler's  rule  of  differentiation; 


5t  +  u 


1  bx±   +  U2  bx2 


+  u.  x~ 

1  ox. 

1 


at 


i  =  1,2  . 


(3O.10 


We  have  also  used  the  summation  convention  in  Eqs .  (3.3.1)  -  (3-3.10  •  This 
convention  will  be  used  throughout  this  chapter.  Furthermore,  the  notation 
of  a  comma  followed  by  a  subscript  i  means  partial  derivatives  of  functions 
with  respect  to  x.. 

The  Eqs.  (3.3.1)  -  (3.3.3),  which  describe  Newton's  principles  of 
motion  for  a  viscous  fluid,  will  be  referred  to  as  time-dependent  viscous 
equations  in  Lagrangian  form.   The  motion  of  the  fluid  thus  described  has  the 
property  of  following  the  position  of  each  single  fluid  particle  for  all  values 
of  t.   However,  there  are  only  four  scalar  equations  in  five  unknowns,  namely, 
p,  U,,  u  ,  p  and  e.   We  need  one  more  equation  to  completely  specify  the 
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problem.   There  exists  no  general  physical  principle  which  would  supply  a 
fifth  equation  to  hold  in  all  cases.   For  perfect  gas,  we  can  use 

p  -  (7  -  l)pe  (30.5) 

where  7  is  the  ratio  of  specific  heats. 

3°  ^°   Equations  in  Streamline  Coordinate  System 

The  solution  of  Eqs„  (3-3«l)  -  (3-3-5)  in  general  describes  a  family 
of  world  lines  in  (x  ,  x  ,  t)-space,  each  line  corresponding  to  one  particle. 
The  projection  onto  the  (x, ,x  )-space  of  the  respective  world  line  is  a 
particle  path.   If,  in  particular,  the  motion  is  stationary,  i.e.,  at  each 
point  of  the  fluid  the  same  thing  happens  at  all  times,  there  can  be  only  one 
family  of  streamlines  in  the  fluid.   The  family  of  streamlines  thus  carries 
most  of  the  flow  information  along  with  it  in  this  case.   It  is  natural  to 
use  the  streamlines  as  one  of  the  fundamental  coordinates  in  the  computation. 

For  stationary  and  psuedo-stationary  flows,  Taub  [16]  introduced  a 
pair  of  independent  variables,  s  and  y.   For  fixed  time  t,  the  curves 
s  =  constant  are  the  loci  of  particles  which  have  crossed  a  given  curve,  say 
y  =  0,  at  times  earlier  than  t.   If  y  represents  the  arc  length  along  the 
curve  s  =  constant,  measured  from  y  =  0,  then  one  may  transform  the  functions 
f  =  f(x  ,x  )  into  functions  of  the  form  f  =  f(s,y)  by  defining 

ax. 

~  =  \.     >  i  =  1*2  ,  (3.^.1) 


and 


s     1 


dx.   u. 

oy  "'  v 


i  =  1,2  (3.^.2) 


where     X.   are  the  components  of  the  tangent  vector  to  y  =  constant; 
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u.  are  the  components  of  the  velocity  vector; 

and 

o 

v  =  u.u.  .  (3-^-3) 

11 

This  defines  the  streamline  coordinates . 

In  addition  to  the  natural  reason  of  using  streamline  coordinates, 
there  are  numerical  advantages .   In  the  point  method  that  we  shall  use  for 
numerical  integration,  a  fluid  particle  can  move  into  any  one  of  four  quadrants 
if  a  rectangular  coordinate  system  is  used.   The  "interpolation"  process  for 
recovering  transport  terms  in  the  equations  of  motion  has  to  handle  all  four 
possibilities  depending  on  the  various  combinations  of  signs  and  magnitudes 
of  the  velocity  components  of  the  particle.   However,  in  streamline  coordinate 
systems,  a  particle  can  only  move  along  a  streamline.   Since  the  streamlines 
are  one  of  our  basic  coordinates,  the  "interpolation"  process  only  need 
consider  two  possibilities.   The  corresponding  "interpolation"  formulas  will 
be  simpler  in  this  case.  A  similar  situation  exists  in  the  cell  method-   The 
amount  of  fluid  moved  out  of  a  given  cell  may  enter  any  one  of  eight  neighbor- 
ing cells  in  rectangular  coordinates.  Using  streamline  coordinates,  the 
outgoing  fluid  from  a  cell  can  only  join  one  of  two  neighboring  cells  along 
the  same  stream  tube.  The  saving  in  computation  is  even  more  apparent  in 
this  case. 

The  second  advantage  is  the  simplification  of  boundary  conditions  in 
the  point  method.   In  streamline  coordinates,  the  streamlines  are  straight 
lines.   Flows  past  inclined  boundaries  are  forced  to  follow  them.   The  stream- 
lines are  parallel  to  the  boundaries.   Values  of  flow  variables  at  a  boundary 

^amline  can  be  obtained  easily.   In  rectangular  coordinates,  the  appearance 
of  irregular  boundary  curves  may  require  complicated  interpolation  formulas 
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to  compute  boundary  conditions.   Similarly,  in  the  cell  method,  if  a  square 
mesh  is  set  up  in  a  two-dimensional  flow  region,  non-square  cells  may  appear 
due  to  the  presence  of  inclined  boundaries.   The  computation  becomes  more 
complicated  there „  Again  the  use  of  streamline  coordinates  can  save  time 
quite  appreciably. 

For  time-dependent  problems,  the  particle  paths  vary  in  time.   They 
can  not  be  used  to  define  coordinates.   The  curves  s  =  constant  defined 
earlier  are  actually  flow  lines  which  now  depend  on  time.   They  shift 
positions  from  time  t  (solid  lines)  to  t  +  At  (dotted  lines)  as  shown  in 
Figure  3.1.   The  correct  interpolation  process  for  flow  variables  at  (£',k') 


y=yo  at  t+At,   P 


s=s  at  t 
o 


FIGURE  3.1   MOVEMENTS  OF  FLOW  LINES  FROM  t  TO  t  +  At 


must  average  the  intermediate  flow  variables  at  (£-l,k),  U,k),  (i-l,k-l)  and 
(i,k-l),  referred  to  as  the  U-point  interpolation .   However,  for  asymptotically 
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stationary  flow,  the  time  dependence  of  flow  lines  approaches  zero.   In  the 
limit,  (i',k!)  coincides  with  (J,k).   The  corresponding  interpolation  process 
at  (-T,k'  )  will  only  involve  the  points  (i-l,k)  and  (£,k),  referred  to  as  the 
2-point  interpolation.   Since  we  are  dealing  with  asymptotically  stationary 
problems,  we  propose  to  ignore  the  movements  of  flow  lines  in  time,  i.e.,  use 
the  2-point  interpolation  in  our  numerical  procedure.   This  simplification 
can  "be  justified  by  the  following  arguments.   First  of  all,  the  asymptotic 
solution  of  our  time-dependent  difference  scheme,  if  one  exists,  is  a  solution 
of  stationary  difference  scheme  and  conversely.   This  is  because  the  flow 
lines  are  stationary  (independent  of  time)  in  the  limit  and  then  the  schemes 
are  identical.   The  stationary  flow  lines  of  the  2-point  interpolation  are 
either  fixed  or  moving  constantly  in  one  direction.   The  latter  is  not 
possible  because  of  the  finite  flow  region.   The  converse  is  true  for  the 
simple  reason  that  stationary  solutions  of  4-point  interpolation  are  solu- 
t ions  of  2-point  interpolation.   Secondly,  the  solution  of  our  difference 
scheme  is  stable  for  small  perturbations  which  implies  that  it  locally 
converges  to  the  solution  of  stationary  scheme  with  respect  to  time.   We  shall 
prove  the  stability  of  the  difference  scheme  with  2-point  interpolation 
process  in  Chapter  IV. 

In  view  of  the  limiting  process  discussed  above,  it  is  possible  for 
us  to  use  the  coordinates  (s,y)  defined  by  the  transformations  (3.^.1)  - 
(3  A. 2)  in  the  proposed  numerical  method.  We  shall  refer  to  them  as  the 
streamline  coordinates,  although  streamlines  do  not  exist  in  time -dependent 
problems  in  the  strictest  sense.   Now,  for  transforming  equations  of  motion 
into  streamline  coordinates,  we  need  the  identity  transformation 

t  =  t  (3.^*0 
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to  take  care  of  the  time -dependent  terms.   Furthermore,  it  is  convenient  to 
introduce  two  new  variables  X   and  OL,    associated  with  X.   by  the  equations  [l6]: 

X2   =  X±X±  ,  (3.4.5) 

and 

tan  a  =  r—   .  (3-^.6) 

The  angle  OL  is   of  course   the   inclination  of  the   curves  y   =  constant   in  the 
(x   ,x   )-plane.      From  [16],   the   following  differential   identities   are  also 
reproduced: 

ui  i  =  if"  chT  '        ^n  not  summed)    »  (3.^-7) 

>  n      y 

df  v  .  df  Uk  ,_   ,    ov 

f    •    =  3T  t"7~   £m  \     -    Y~   £M     77"      >  (3-4.8) 

,J       oy  Un     jk  k       ds     jk  Un 

cp  =   e..u.     •    =;f  (I1"  ^~)      >  (3.^.9) 

u.f    .=y|     ,  (3.i+.10) 

and 

X.f    .    =^f     ,  (3.4.11) 

where  U.    =  u.\.      ,  (3-4.12) 

til 

is   the  velocity  tangential  to  the   curves  y  =  constant,    scaled  by  \; 

U     =u.€.A.      ,  (3A.13) 

n  1   ij   J 

is  the  velocity  normal  to  the  curves  y  =  constant,  scaled  by  X, 
Snd  £11  =  €22  =  ^  '£12  =  "*21  =  1   ' 


a  ^ 
COS    p   =r— 

ds 

> 

da 
Sy  = 

sin   P  d£ 

a.    as 

J 
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If  the  angle  (3  is  defined  as  the  angle  between  the  normal  to  the 
curve  y  =  constant  and  the  tangent  to  the  curve  s  =  constant,  then 

U  =  \v  sin  P,  (3.4. l4) 

and 

U  =  A.v  cos  p.  (3.4.15) 

The  angle  P  is  taken  to  he  positive  when  the  direction  from  the  normal  to 
y  =  constant  to  the  tangent  to  s  =  constant  is  counter-clockwise.  It  can 
also  be  shown  that 

(3.^.16) 

(3.4.17) 

=  v  cos (a  +  p  -  — )  =  v  sin  £  (3.4.18) 

and  up  =  v  sin(a  +  p  -  — )  =  -v  cos£  (3.4.19) 

where     %   =  a  +   P  .  (3.4.20) 

Thus  the  angle  (a  +  P  -  — )  gives  the  inclination  of  the  curves  s  =  constant  in 
the  (x  ,x  ) -plane. 

Due  to  the  presence  of  viscosity  terms,  we  shall  need: 

(ui,Aj  ■  (\,k>,i +  Ki\,fu\i  -  K*\i +  ^u\t  (3-u-21) 

Using  the  Eqs.  (3.4.1)  -  (3.4.21),  we  can  transform  the  Eqs.  (3-3.1) 
■3)  into  the  new  coordinate  system  as  follows: 

(l)  Continuity  Equation 

dU 
11+^^  =  0  ,      (n  not  summed)  (3-4.22) 

n 


u 
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or 


p  =  -pi  , 


(5.^.25) 


where 


du 

v    n 

U  dy  ' 

n  " 


(5-^.2^) 


(2)  Momentum  Equations 


"h^KiK^hK^A'^A 


i  =  1,2  . 


(5^.25) 


Multiplying  Eqs.  (3.^.25)  "by  u.  and  summing,  we  have 

pu .  u .  =  u.u .  (cpe  .  . )  .  +  —  u.u .  (un  .  )  .  -  u .  p  .   , 
K  1  1   ^  x       xljl        3   x  k,k  ,1    1  ,1 

where  we  have  used  Eq„    (3.4-21).      Using  Eqs.    (3«^<-7)    _    (3-^-H)   and  simplifying, 

the  above  equation  becomes 

2n_      ,vU.    v  1  dU 


Or, 


Pw    =   u.(—  ^   -    -77-  k^-J    +    —  [IV   ^—    (77-  xT~)    -    V    AT 

^VU     ds         U     dy         3         dy     U     dy  dy 

n  n     J  '  J         n     ^ 


v     ,dcp  t   dou         ^        dt        di 


U      vds         v 

n 


5  »  *y 


(5.^26) 


Multiplying  Eqs.    (3 A. 25)   by  \.    and  summing,   we  have 


pX.u.    =  ^(cpc.p^   +  -  ^(u^  .    -   ^.    . 


Or, 


4       d 


v      Ut   dcp         2  dcpN 
p\(|v  cos   p  +  v  sin  p)   =  n  ~  (-  g  -  \     ^)    +  -  u  ^  ( 

n 


dU  . 

y     n\  dp 

3  ^  ds  VU    oy~;  '  ds   ' 
n 


Finally,   using  Eqs.    (3.K°lk)    -    (3.^.20)   and  Eq.    (3.J+.26),   we   can  write  the 
above  equation   in  the   form 


h  .     rdt       jt   d\J/\  Jl  M  ..  ^£       _i  ^P    • 


PV    =  5  ^    <£  -  T  d7}    "   ^ 


v    dy       ds         v    dy 


(3.^.27) 
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(3)   Energy  Equation 


>(e  +  5  Yl 


-    u(u.      .U.       .     +    U.      .U.      .      -    —  U .      .LL  ) 


u .(u.     . )    .   +  u  .  (u  .    .  ) 
J      1,1    , 3  J      J,i    ,ij 


-  —  uu .  (u.    ,  )    . 
3        i     k,k  ,i 


-  pu. 


1,1  1,1 


u.p    .    -   p(e  +  -  v   )u.    . 
,i  2  j, j 


Or, 


(pe)      +  -   (pv    ) 


k     .2  2  ,  y^6     dv  d£     OVs 


+   H 


v   xr  +  u.(\|r    .   +   (cpe.J    J      -  -  uv  3- 


p^   _   v  ^  _   p(e  +  -  v   )\|r    . 


Or, 


(pe)'    = 


4    .2  2        ■       ,d£   dv        d£    dvx        If        ot       v      ,dcp 


u, 


.3 


v   oy;_ 


-  v  J*  -  p\|r  -    (-  pv   )      -p(e-+  -  v  )t   . 

Substituting  Eqs .    (3.^.23)   and   (3.^.26)    into  the  above  equation  and  simplifying, 
we  obtain 


pe  =  u 


k-    . 2  2        ,     ,d£   dv        as   dvx 


-  pV 


(3.^.28) 


Let  us  define 


and 


U 

V 


(3.^.29) 
(3.^.30) 


We  then  have  from  Eq.    (3.I+.2I+) 
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,        1  a      /,  Q\        v  hr\        av 

•ty    _    _   __    t\v    Cos     PJ     =    —   -r-1   +    t— 

t\  oy  t)  ay       ay 


where  we  have  used  Eqs .  (3  A.  14)  -  (3-4.17)  and  Eq.  (3.4.20) 


(3.4.31) 


Similarly ,   we  can  write  cp  in  terms  of  r\   and  £  as  follows 


1  /dv     a£    «.  cVv\ 


The  differential  equations  for  t]  and  £  are  simply, 


ay  "  as  "  *•  dy  ' 


(3.4.32) 


(3.4.33) 


and 


oy  "  n  oy 


(3.4.34) 


Substituting  Eqs „  (3.4.33)  -  (3*4.34)  into  Eqs.  (3.4-31)  -  (3.4-32)  respectively, 
we  get 


♦  -?&-■«  !>♦!' 


1     /OV  u,    GV-, 

*  =  n  (^  "  s  ¥} 


as 


Finally,  the  equations  of  motion  in  terms  of  r\   and  £  are: 


p  =  -pt   , 


pv  =  ^   (- 


acp\     4 ,,  M    ap 


oy         3       ay       oy 


pvl  -  3  ti  ^as     Q  ay;     *-  ay 


n  ay 


1  ap 

T)    3s~     ' 


and 


=  n(7-l) 


4  .2      2     .    ,ai  av     as  av^ 
3  *  +  *   +  Uv  (¥  as" "  ai  ¥} 


7pt      , 


where  we  have  used  Eq.  (3.3=5)  to  obtain  Eq.  (3»4.4o) 


(3.4.35) 
(3.4.36) 

(3.4-37) 
(3.4.38) 

(3.4.39) 
(3.4.40) 


Eqs.  (3.4.33)  -  (3.4.40)  form  a  set  of  eight  first  order  partial  differential 
equations  in  eight  unknowns  p,  v,  £,   p,  t\,    t>,   cp  and  ^o   They  describe  the 


time -dependent,  two-dimensional  viscous  fluid  flow  phenomena  in  streamline 
coordinates  when  appropriate  boundary  conditions  are  specified. 


CHAPTER  IV 
NUMERICAL  SOLUTION  OF  VISCOUS  SHOCK  REFLECTION  PROBLEMS 

4.1.   Stability  of  One -Dimensional  Difference  Equations 

In  this  section,  we  are  only  concerned  with  the  stability  of  one- 
dimensional  equations ,   We  thus  have 

Un  =  v,      Ut  =  0.  (4.1.1) 

This  leads  to  the  conditions  that 

r\   =   1,      £  =  0  ,  (4.1-2) 

by  definition.   Substituting  Eqs .  (4.1.1)  -  (4.1.2)  into  Eqs.  (3.4.31)  - 
(3°4c38),  we  get  the  following  one -dimensional  viscous  equations  of  motion 

P  =  -pf  >  (4.1.3) 

-  -  M  - 1  >  <^-« 

4        2 
p  =  -  n(7-l)t  -  7pt  ,  (4.1.5) 

and 

*■¥•  (^6) 

For  the  sake  of  stability  analysis,  it  is  convenient  to  eliminate  \|/  from  the 
above  equations  by  using  Eq.  (4.1.6).   We  then  have 

P=-p|,  (4.1-7) 

pv^.g§-|,  ft.1.8) 

ana  p  = ; .  r-ixfr)   -  ^p^  •  (k-1'9) 

For  evenly  spaced,  mesh  points,  we  employ  central  differences  to  approximate 
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space  derivatives.   For  the  time  derivatives,  we  use  the  one-sided  difference 
approximation.   Writing  Eqs.  (4.1-7)  -  (4.1-9)  in  difference  form,  we  have 

n+1    n        n      n 

■  p»     n  Vl  ~   Vl 


p.  ^~ >  (4.i.ioj 


At        HH  2Ay 

n+1    n        n     0  n    n      n      n 
n  Vl    -  Vl    4   Vl  ~   2Yt    ~   Vl   Pl+1  -  Pl-1  ,.-_.* 

Pi   At"  =  5  ^  ^2~ 2Ay" '  (^1'11) 

n+1  n  n  n        2  n  n 

Pl       '  Pl       4     .     lN/i+l  "  VI-1N         „  n  vi+rVi  ,,  » 

and  ______  ,_M.(7_1)( ^y        )      -  7Vi  ^ (  +  .1.12) 

Assuming  that  a  solution  of  Eqs.  (4.1.10)  -  (4.1.12)  can  be  expressed 
in  terms  of  Fourier  series,  let  us  examine  the  conditions  under  which  a 
typical  term  from  the  series  can  be  a  solution.   Take  as  the  representative 
terms  respectively 

n       _   ikJ  n  n    .  __x 

Pjl    =   p  +  Spe   r   ,  (4.1.13) 

n       _   iki  n  ,,    _  _ ,  v 

v  =  v  +  Sve   r   ,  (4.1.14) 

n       _   ikl  n  , ,  n  ,  _ N 

and      p  =  p  +  ope   r      ,  (4.1.15) 

where  p,  v  and  p  are  the  true  solution  of  Eqs.  (4.1-7)  -  (4.1-9) • 

Substituting  them  into  Eqs.  (4.1.10)  -  (4.1.12)  and  simplifying,  we  obtain 

(r-l)Sp  =  -ap(el  -  e  X  )6v  -  (^  At)6p  ,  (4.1.16) 

p(r-l)8v  +  (vAt)Bp  =  I  in(eX   -2+e"1  )&v-a(— |-S — -)op    (4.1.17) 

and      (r-l)Bp  =  |  n(7-l)o  *   (£_^|__)&v  -  7pa(£_^|-_)5v  -  7  ^  At6p   , 

(4.1.18) 
ire    0  ==  p^  ,  (4.1.19) 


and 


T  = 


At 

Ay2 


(4.1.20) 


Assuming  \x,  At  and  Ay  are  small  and  of  the  same  order  of  magnitude,  we  can 
hold  a  and  [IT  fixed  and  let  At  and  [ia  -»  0  in  Eqs.  (4.1.16)  -  (4.1.18).  In 
the  limit  we  get  the  following  system  of  linear  algebraic  equations  in  5p, 
&v  and  op 


(r-l)bp  +  (lap  sin  k;bv  =  0  , 
8 


p(r-l)  +  ~  u:  1  ■  cos  k) 


5 v  -»-   i  a  s  1  n  k  I  op  -  C  , 


and 


L7po  sin  k)8v  +-  (r-l'Sp  =  0 


(4.1.21) 
(^. 1.22) 

(4  1.23) 


A  necessary  and  sufficient  condition  that  Eqs   (4.1  2L   -  (4.1.23)  have  a 
non-trivial  solution  for  5p.,  ov  and  op  is  that 


r-1 


iap  sin  k 


C 


0      I  r-1,)  f  —  |J.t{1  -  cos  k'     i.a  sin  k 


i?pa  sin  k 


(r-1 


=  0 


(4:1. 2*0 


Or 


1   1  >2   8 

J r-1)   +  -  q  - 

5   P 


5  2   2 
L  -  ccs  k)(r-lj  ■•-  c  a  sin  k 


(4.1.25! 


where 


P 


I U.1,26] 


In  obtaining  Eq<  (4.1.25),  we  have  discarded  tne  asymptotically  stable  root 
r  =  I.   The  solution  of  Eq.  (4.1  25'  can  be  written  in  the  form 


r   =  1  -  nil 


I  g  2   2   2  ' 

cos  k)  -  A/a  (l-cosk)  -p  sin  k  , 


wl  e  re    a. 

and      p. 


4   1 

3  "  e 


ca 


(4.1.2?) 
(4.1.28) 
(4.1.29) 
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In  order  that  solution  (4.1.13)  -  (4.1.15)  does  not  grow  with  time, 
it  is  necessary  that  |r|  <  1.   This  condition  imposes  restrictions  upon  the 
parameters  At  and  Ay.  We  will  discuss  the  stability  conditions  in  two  cases: 

CASE  I.   a2(l-cosk)2  <  p2sin2k 

Then  r  is  complex  and 

|r|2  =  1  -  2a(l-cosk)  +  p2sin2k 

=  1  -  (l-cosk)[2a-  P2(l  +  cos  k)]  .  (4.1-30) 

Thus,  in  order  that  |r|   <  1,  it  is  required  that 

2a   -  p2(l  +  cos  k)  >  0  . 
Using  Eqs.  (4.1.19),  (4.1.28)  and  (4.1.29),  we  finally  get 

4 

-   7P 

where  we  have  replaced  (l+  cos  k)  by  its  maximum  value  2. 

Case  II.   a2(l  -  cos  k)2  >  p2sin2k 


At  <  -* (4.1.31) 


In  this  case,  r  is  real  and  -vj-  is  proportional  to  sin  k  so  that 


dr 
dk 

|r|  attains  its  maximum  value  at  k  =  0,  it.   For  k  =  0,  r  =  1.   No  restriction 
is  imposed  on  At  and  Ay.   For  k  =  jt,   we  have 

r  =  1  -  4a  .  (4.1.32) 

For  stability,  we  require 

a<|  .  (^.1.33) 

This   gives 

At  <  td^dl  (.k.i.ik) 

5» 
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For  sufficiently  small  Ay,  condition  (4.1„3^)  dominates.   We  will  use  it  to 
determine  the  time  step  for  numerical  integration, 

k . 2 ,   Nume r i cal  Solution  of _One -Dimensional  Viscous  Shock  Problem 

Lmagine  a  long  pipe  with  one  end  closed,   A  uniform  fluid  flow  is 
"being  supplied  continuously  from  the  other  end.   When  the  fluid  particles 
meet  the  closed  end,  disturbances  are  created.   A  shock  front  is  generated 
which  moves  away  from  the  closed  end  as  shown  in  Figure  k<l 


FIGURE  k.l        ONE -DIMENSIONAL  VISCOUS  SHOCK  PROBLEM 


In  the  absence  of  dissipation,  the  density  p  and  the  fluid  velocity  v  are, 

at  a  given  instant,  as  shown  by  the  solid  curves  in  the  graphs,  whereas  in  the 


-1*0- 

presence  of  dissipation,  they  are  as  shown  by  the  dashed  curves.  The  shock 
front  is  smeared  out  when  viscosity  effect  is  included.  Analytically  we 
can  obtain  the  solution  in  closed  form  if  the  heat  conduction  terms  are 
ignored  [ 17] • 

For  the  numerical  solution  of  this  problem,  we  apply  the  finite 
difference  approximation  to  the  Eqs .  (4.1. 3)  -  (4.1.6)  for  L  equally- spaced 
mesh  points  along  y-axis  and  get 


-n+1 


n   AJ_  /  n  ,n\ 


I 


1  Tl 


(4.2.1) 


,n              n 

-,           1+-          1 — 

-n+1 

n       At 

4               2               2 

Vl 

z  v     +  ■ — 
1          n 

3^              Ay 

1+1 


'1-1 


2Ay 


(+.2.2) 


and 


—n+1    n 

Pfl   =  Pf  +  At 


nN2 


n ,  n 


^(7-D(t;r  -  rP;> 


(4.2.3) 


where 


n 

<4 


n      n 
v    -  v 
1+1    1 


Ay 


(4.2.4) 


and 


n    n 
v  -  v 
1    1-1 


Ay 


(4.2.5) 


n   1  /  ,n 


, n   1  / ,n      n   N 
*!  -2  (*  1+  *  1} 


(4.2.6) 


!+—    !• 
2      2 

The  reason  for  using  half  mesh  values  of  \|f  is  to  keep  uniform  error  in 

difference  equations. 

Since  p    ,  v   ,  and  p  '  are  the  intermediate  flow  variables  of 

the  1th  mesh  point  at  ( n+1 ) -time  cycle,  we  need  to  combine  them  with  the 

of  appropriate  neighboring  mesh  point  for  the  "interpolation"  process. 
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The  weighted  linear   interpolation   is  performed  according  to  the  following 
formulas : 

CASE  I.        v^  >  0 

n+1   .— n+1    — n+1  ,,     x 

Pi   =  Api-i  +  Bpi   '  ^-2-T) 

n+1   A— n+1   -n— n+1  ,,  „  ON 

Y*    =  AVi  +  Bvi    -  ^-2-8> 

and  p^+1  .  Ai£  +  Bp^+1  ,  (4. 2. 9) 

v  At 

where         A   = — — —  ,  (  +  .2.10) 

n         n 
Ay  -  v^_x  At  +  v^  At 

n 
Ay  -  v    At 

and  B  =  ■ =^±- .  (4.2.11) 

Ay  -  v^_x  At  +  v'J  At 

For  allowing  mesh  points  to  move  in  opposite  direction,  we  need 

CASE  II.   v^  <  0 

n+1    — n+1   — n+1  /,  _  __\ 

pi    =A^i    +Bpi+1  '  (^2-12) 

n+1   .  —n+1   t,— n+1  / ,  _..,-,  \ 

and  p^+1  =  A^+1  +  B?£  ,  (4.2.14) 

Ay  +  v"   At 

where         A  =  — ,  (4.2.15) 

Ay  +  v^+1  At  -  v^  At 

and  B  -  - —   .  (4.2.16) 

Ay  +  v°w  At  -  v«  At 

The  reasons  for  choosing  above  formulas  are  quite  clear  geometrically  as  shown 
in  Figure  4.2. 
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^V-VJL,At       JC?M^ 


v£>0 


-v^n+1<o 


FIGURE   4.2        INTERPOLATION  AT  MESH  POINTS 


(n  +  l)-TIME  CYCLE 


(n)-TIME  CYCLE 


in  the  stability  analysis  of  Section  4.1,  we  have  only  shown  that 
the  difference  scheme  is  stable  in  Lagrangian  form.   We  need  further  proof 
that  it  is  also  stable  after  the  interpolation  process.   Let  us  consider 
Case  I  for  simplicity.   Let  f  be  any  flow  variable  being  interpolated.   From 
Eqs.  (4.2.7)  -  (4-2.9),  we  can  write 


„n+l 


*?" 


^i+l 


+ 


^  '  Vl  At     ?n+l 
T    ""    n  ..  i 


Ay  -  v*   At  +  v*  At  "i"1         "_i 


'A-l 


Assume  v»  =  v^  -  v.  And  using  Eqs.  (4.2.1)  -  (4.2. 5),  we  have 


-n+1 


and 


f^+1  .  fj  +  At(f)* 


(4.2.17) 


(4.2.18) 


(4.2.19) 
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Eq.  (4.2.17)  now  becomes 

Substituting  f  in  the  form  of  Eqs .  (4.1.13)  -  (4.1.15)  into  the  above 
equation  and  simplifying,  we  have 

r  _  2**  (e"lk  +    At     fn   )  +  Ay  -  vAt  (1  +    At     fn) 
Ay  K*  +  R^  ik*  n  Lt-lJ  Ay     l     ^  ikl  n  'l^ 

J  5fe   r  J         5fe   r 

fn 

»!_£*(!_  e-lk  -  ^   -.;.  ■  )  +  0(At2)  .  (4.2.21) 

Ay   v  v  c_  lki  ny     v    '  K  ' 

ofe   r 

In  the  limit  as  At  -»  0  for  fixed  small  Ay,  r  approaches  the  asymptotically 
stable  root  r  =  1.   For  all  small  At,  |r|  <  1.   Thus  the  interpolation  process 
is  stable  locally.   This  analysis  together  with  At  satisfying  (4.1.34)  for 
sufficiently  small  Ay  give  the  stability  of  difference  scheme  with  interpola- 
tion in  streamline  coordinates. 

The  boundary  conditions  at  the  open  end  are  simply 

nn  nn  nn  f\.   n  nn\ 

P2  =  P-l      ,  v2  =  vl      '    p2  =  Pl  (4.2.22) 

to  ensure  the  continuity  of  the  fluid  flow.  At  the  closed  end,  we  use 

nn  nn  nn  /\.   n  ni\ 

pL  =  PL-1    '     VL  =  "VL-1    '     PL  =  PL-1   '  (^.2.23) 

The  initial  flow  values  at  all  mesh  points  are  uniform.   They  are  taken  to  be 

p°  =  1.0    ,    for    I   =  2,3,..., L-l   , 

v°  =  1.0    ,    for    1  =  2,3,..., L-l   ,  (4.2.24) 

and       p°  =  4.0    ,    for    I  =  2, 3,..., L-l   . 

The  value  of  p  is  chosen  to  give  supersonic  incoming  flow.  Other  constants 
required  for  numerical  integration  are 
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L  =  50,   |  [i  =   0.2,   7  =  1.4  .  (4.2.25) 

3 

And  At  is  chosen  to  satisfy  (4.1.34). 

On  Figure  4-3,  we  plot  the  results  obtained  with  different  combina- 
tions of  At  and  Ay  all  satisfying  the  stability  requirement  (4.1.34).   The 
wavy  pattern  of  the  velocity  and  density  behind  the  shock  front  observed  by 
Richtmeyer  [12]  and  Crocco  [5]  shows  quite  clearly  in  the  graphs  in  cases  B 
and  C.   It  was  shown  analytically  [5]  that  these  waves  appear  whenever  the 
space  interval  becomes  larger  than  half  the  shock  thickness.  A  simple 

explanation  is  that  the  inaccuracy  of  the  finite  difference  approximation 

p 

0(Ay  )  comes  in  when  Ay  is  too  large. 

4.3'   Stability  of  Two-Dimensional  Difference  Equations 

The  stability  analysis  in  this  section  is  similar  to  the  one  used 
in  Section  4.1.   Here  we  have  to  use  the  basic  equations  (3-4-37/  -  (3A.40) 
in  their  full  generality,  however.   Those  who  are  only  interested  in  the 
resulting  stability  condition  may  skip  the  detailed  derivations  to  page  50. 

Let  us  impose  a  rectangular  mesh  upon  the  two-dimensional  region 
under  consideration  and  use  the  same  difference  approximations  as  in  Section 
4.1.   It  is  again  convenient  for  stability  analysis  reasons  to  eliminate  cp 
and  \|f  from  the  Eqs .  (3-4.37)  -  (3-4.40).  We  then  take  the  representative 
terms  from  the  Fourier  series  expansion  of  each  flow  variable  and  substitute 
them  into  the  resulting  difference  equations.   To  simplify  the  resulting 
equations,  we  neglect  terms  of  higher  orders  in  differentials  of  flow 
variables.  After  some  tedious  manipulations  and  simplifications,  we  get 

T6p  =  .  £Z  (S  .  £Y)B£  .  pY6v  (4-3-1) 
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FIGURE   k.J>        WAVY   PATTERN  BEHIND  THE  SHOCK  FRONT 
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pTSv  =  -^   (S^Sv   -  VTiSYbl    -    £SY8v)    -  -^  £(YSSv   -  vt^&I    -^Bv) 

T)  Tl 

+  uj-   (YS   -    £?)&£    +  ?6v]    -  YSp  (4.3-2) 

pvT&^    =  -^   [^  (S2   -    £SY)B£   +   SY   5v]    -  -^  £[Z   (YS    -£?)&£    +  Y^Bv] 

-  ^   (YSBv   -   vtJY25£    -    ^Bv)    +  -^  YSp   -  -  SBp  (^.3.3) 

and  TBp   =    (7-l)[2nQ\|f(^   (S-£Y)B|   +  YSv)    +  ^   (SBv   -   vi]YS£    -    £Y8v) 

+   kavii   SBv  +  v  YS£    -   g  YBv   -  v  SB|) ] 
y  s  s  y 

-  ?P[^   (S   -^Y)8^   +   YBv]  (4.3-4) 

where  ij.     =     —  p.   ,  (4.3-5) 

T     =     (r  -   1)/At   ,  (4.3-6) 

1  Sln   i  (4.3-7) 

(4.3.8) 

(4.3-9) 
(4.3-10) 

(4.3.11) 

Ah 

Ah  =  As  =  Ay  for  simplicity,  (4.3.12) 

and      a  subscript  s  or  y  means  partial  derivative  with  respect 
to  s  or  y  respectively. 

Substituting  Eqs .  (4.3-6)  -  (4.3.12)  into  Eqs .  (4-3.1)  -  (4.3-4)  and  collecting 
terms,  we  have 

(r  -  l)8p  +  ^2.  i(sin  l    -  £  sin  k)&£  +  pa(i  sin  k)&v  =  0, 

(4.3.13) 


0   - 

Ah     ' 

Y   = 

i  sin  k 
Ah 

s2  = 

2(cos£  -  1) 

2     ' 
Ah 

Y2  = 

2(cos  k  -  l) 

2      ' 
Ah 

RY  - 

^rr,     sin  k  sin  I 

and 


-fcT- 

i(r   -   l)Bv  +  ^r  x  i[2(cos    £    -   l)5v  +  vt]   sin  k  sin   £5| . 

+   £   sin  k  sin   Hbv]    -    £[sin  k  sin   Ibv  +  2vr](cos   k  -   l)6£ 

+   2£(cos   k  -   l)5v]  ?•  -  \i  t[2(cos   k  -  l)Bv   -  —  (sin  k  sin   I 

+  2£(cos  k  -  l))5£]   +  a(i  sin  k)sp  =  0   ,  (k.J.lk) 


pv(r  -  l)S|    -  -  1u  [-  (2(cos    I   -  1)   +   5   sin  k  sin   l)bi 


v 


-  sin  k  sin   Ibv]   +  p.   £[2(cos   k  -  l)Sv  -  —  (sin  k  sin   I 

+  2£(cos   k  -   l)8|]    -  u[sin  k  sin   i&v  +  2vt](cos   k  -  l)6£ 

+  2£(cos   k  -   l)Sv]f  -  •£-  ai  sin  k6p  +  -  i  sin  *6p  =  0,  (k.J.lj) 

(r   -  l)6.p  -   eri  -j(7    -   l)[2u  \|r(-  (sin   I   -    £   sin  k)5|  +   sin  k6v) 

+  2  «-  (sin   Ibv  -  vr\   sin  kS£    -    £   sin  k5v)   +   kp.v   (£      sin   £5v 

+  v     sin  k&£    -   |      sin  kSv   -  v     sin   i&Ol 
s  s  y 

+  7p[-  (sin   I    -   £   sin  k)&|   +  sin  k&v] f     =  0.  (4.3-16) 

Eqs.  (4.3-13)  -  (4.3-16)  form  a  system  of  linear  equations  in  5p,  6v,  6£ , 
and  5p.  A  necessary  and  sufficient  condition  that  they  have  a  nontrivial 
solution  is  that  the  determinant  of  the  coefficient  matrix  A  -  (a .  . )  be  equal 
to  zero .   The  elements  of  A  are i 

an  =  r  -  1  , 

a,   =  pai  sin  k  , 

an _  =  ^—  ai(sin  I    -    £  sin  k)  , 

13     T]  v  3         '  ' 

aik    =  0  ' 
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a21  -  °    ' 


=  p(r       l)    -  \  i2u.[(cos    £    -  l)   +  |  sin  k  sin   I] 


'        „2 


a23 


-  u.£[2£(cos   k  -   l)   +   sin  k  sin   I]    +  2(jl  t}    (cos   k  -   l)  f  , 

=  ■ —  -m   [sin  k  sin   I   +  2£(cos   k  -   l)]    -   u.   sin  k  sin   I 

-  2u.£(cos   k  -  l)Y  , 


a2k  =  0±   sin  k   ' 


a51  =  0  , 


T 

a3£  =  n 


a53 


■j   u.     sin  k  sin   I   +  2(i  £(cos   k  -   l) 
-  u.[sin  k  sin   I   +  2£(cos   k  -   l)]f  , 
=  pv(r   -   l)    -  —  iu.   [2(cos    I   -  l)  +   £   sin  k  sin   £] 

+   u.   £[sin  k  sin   I  +   2£(cos   k  -   l)]    +  2u.t]    (cos   k  -   l)  f  , 


a 


=  —  (sin   £    -    £   sin  k)    , 


3^  '     n 

V  =  °  ' 

a,      =  7pai  sin  k  -  (7   -  l)ai[2u.  tysin  k  +  -"•  (sin   £   -   £sin  k) 

+   Vv(£      sin  I   -  i      sin  k)  ]    , 

a,      =  --21  ai(sin   i  -   £   sin  k)    -    (7   -   l)ai[2|_i  ij/  —  (sj_n   1    _   £sin  k) 

-   2p.cpv  sin  k  +   4u.v(v     sink-v     sin   i)]    , 

s         y 

%k   =  *  -  1  • 

Discarding  the  asymptotically  stable  roots  r  =  1  as  before  and 
assuming  0  -»  0  faster  than  t  for  sufficiently  small  Ah,  we  are  left  with  the 
equation 
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(r   -   l)2    -   2  ^~  1  [(cos    k   -   l)(£2+q2)    +    (cos    I    -   1 

P1! 


2   2 


+   £   sin  k  sin   i)]    (r-l)    +        ^, 

P   n 


16    /  i  A2/*-2  2\ 

—  (cos   k   -   1;    (£      +  T)    ) 


+   4(cos   k  -   l)(cos    I    -  l)   +    ^  sin  k  sin   l)(-  ^     +q11    ) 


16  2        n  ' 

+  —  (cos    I    -   1  +  ^sin  k  sin   t)      -  4r-  (2(cos   k  -   l)£  +   sin  k  sin   i)' 


4   /  ,         ,v2    .2   2 

+   -   (cos    k   -    1)       I    T] 


=     0 


(4.3-17) 


The  roots  of  Eq.  (4. 3° 17)  can  be  shown  to  he 


■Y>  —  _  -  ~ 

A1      2  2 

Ah     prj 


J  B  T  i  Vb2  +   C 
.3  5 


(4.3-18) 


2  2> 

where  B  =   (l   -   cos   k)(£      +  T)   )  +   (l  -   cos    i    -    £   sin  k  sin   £)    ,  (4.3*19) 


2  2  2  2 

and  C  =  T]      sin  k  sin   I    -   4tj    (cos   k  -   l)(cos    I    -   l) 


(4.3.20) 


For  stability,   we   require 


At     < 


(Ah)2pi]2 


f  b  ;  i  ^B2  +  c 


(4.3.21) 


But 


,^2„2k.2i/        2k  2   1        ,^      ,     „ 

C   =  lor)   sm     —  sin     —  (cos     —  cos     —  -  1)      <     0 


(4o3-22) 


Therefore,   we  can  simplify   (4»3°2l)   and  get 


At      <        WW        < 
^  |i  Max(B) 


f  Bf  \  Vb2  +  C 


(4.3.23) 


The  extrema  of  B  occur  at  k  =  0>*  and  i  =  0,it.   We  have  the  following  cases  : 
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(1)  k  =  0  ,  I  =  0  , 

(2)  k  =  0  ,  i  =  «  , 

(3)  k  =  Tt  ,  I  =  0  , 

(k)     k  =  jt  ,  £  =  a  , 

Condition  (4.3-23)  now  becomes 


B  =  0  , 
B  =  2  , 


(minimum) 


B  =  2(£2  +  T!2)  , 

2    2 
B  =  2(£  +q   +  l)   (maximum) 


At  < 


/      (Ah)2pn2 
-  15  772    2    _v 


(h-3.2k) 


Since  we  neglected  the  term  C  in  deriving  (4.3*23) >  the  stability  condition 
(4-3.24)  may  be  a  little  too  restrictive.   However,  it  does  give  us  a  reasonable 
upper  bound  for  At  in  terms  of  Ah.   Comparing  to  the  corresponding  one- 

dimensional  stability  condition  (4.1-34),  we  can  get  £  =  0,  T]  =  1  and  £  s  0. 

2    2 
The  maximum  of  B  now  takes  on  the  value  2(£   +  t]  )  =2.   (4.3-23)  then 

becomes 


At  < 


p(Ah)' 


-   16 


(4-3.25) 


I" 


This  gives  a  smaller  upper  bound  than  (4.1.34)  because  of  the  simplification 
made  earlier.   For  the  proof  of  stability  of  interpolation  process  in  this 
case,  it  is  identical  to  that  given  in  Section  4.2  since  the  interpolation 
formulas  are  the  same.   We  therefore  have  the  stability  of  the  two-dimensional 
difference  scheme  with  interpolation. 


4.4.   Numerical  Solution  of  Regular  and  Mach  Reflection^ 

In  this  section,  we  shall  describe  a  numerical  procedure  of  solving 
shock  reflection  problems.   To  simplify  the  computation,  we  consider  a  two- 
dimensional  rectangular  box  (Fig.  4.4)  where  uniform  fluid  flow  comes  in  from 
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the  left  and  exits  at  the  right.   The  top  and  bottom  boundaries  are  assumed 
rigid  and  smooth.   The  boundary  layer  effects  will  be  ignored  completely. 
The  upper  wall  can  be  thought  of  as  the  dividing  streamline  of  two  identical 
shocks  incident  on  each  other  as  explained  in  Section  2o2. 


UPPER  WALL 

///////////////////////////////////////////// 


////// 
LOWER    WALL 

FIGURE  k.k       TWO-DIMENSIONAL  SHOCK  REFLECTION  CONFIGURATION 


One  way  of  generating  a  steady  incident  shock  in  the  flow  is  to 
place  a  straight  half  wedge  with  angle  go  on  the  bottom  boundary  of  the  box. 
Flow  will  then  be  disturbed  starting  at  the  tip  of  the  wedge.   If  the 
incoming  flow  is  supersonic  and  tne  wedge  angle  is  small,  a  fairly  straight 
attached  shock  will  appear  as  shown  in  Figure  h.K.      Using  this  simple 
arrangement,  desired  incident  shock  angle  a     and  its  strength  ^  can  easily 
be  achieved  by  fixing  the  half  wedge  angle  u  and  the  incoming  flow  Mach 
number  MT. 


-52- 

In  the  present  calculations ,   a  square  mesh  is  set  up  in  the 
rectangular  box  in  streamline  coordinates.   The  finite  difference  equations 
applied  at  each  mesh  point,  derivable  from  Eqs .  (3 •^•37)  ~  (3-^-^0),  can  be 
written  in  the  following  form 


n+1    n        n    n 


(k.h.l) 


9 


9 


n+1    n      At 

v„    =  v.   +  

l,m  I  ,m         n 


l,nn 


,m-g 


n, 


Ah 


-  5 


9,  i  _cp,  i 

n    *+2'm    £-2-m 
!,m  "      Ah 


£,m   £,m 


+  M> 


,n        n 

.  1      .  1  n  n 

i+— ,m    ^-7T>m  P.  ,    -  P„  -, 

2' 2*  l+l,m    l-l,m 


Ah 


2Ah 


(4. +.2) 


, n+1   , n 

I.      =  £.   + 


At 


£,m    i.m    n    n 
P  -   v 


I-1, 


Tl 


n       ,n 
+    1"  +    1 

k"      Ah 


£,m  £  ,m   l  ,m 


9  !   -cp 


,n       ,n 
*  1   "  +  1 

^£,m        Ah      '  "  ^  ~      Ah 


.n    n        n  n        n 

^l,m  Pl+l,m  "  Pl-l,m     1   Pl,m+1  "  pl,m-l 
n         2Ah       "  n         2Ah 


'i,m 


'i,m 


(^•3) 


and 


n+1 
~*  l,m 


C  +  *  fa-v  [j  <*?fB)2  +  KJ 


»"  *  n  "  ii 

liv11        /    *+l,m     l-l,m  Vl,m+l~Vl,m-l 
Vi,m  2Ah  2Ah 


»n  n  n  n 

/,m+l~Sl,m-l  vf+l,m"V*-l,m 
2Ah  2Ah 


n         in 
P£,m     £,mj    . 


(+.+.+) 
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In  order  to  keep  uniform  error  in  the  above  difference  equations,  we  use  half 
mesh  values  of  cp  and  \|/  in  the  corresponding  difference  quotients  „  These  half 
mesh    values   can  be   computed  from  Eqs .    (3. 4. 35)    -    (3 <^«36)   as   follows 


t 


n 
i+-,m 


2> 


m 


n 


i 


l+-,m 


£+^,m 


-y .  i 


n  n 

v  -v 

£+l,m     I ,m 

Ah 


(^•5) 


and 


f 


£  ,m+ 


1 


£  ,m+- 


£,nn- 


£,m+2 


<oD 


i,m+ 


i    3y 


,m+~ 


n  n 

v  -  v 

£,m+l  £,m 

Ah 


(h.h.6) 


where  v 


n 
£+2>m 


1    /   n 


(v 

V       1  +  1, 


m  £  ,m 


(+.^•7) 


V 


£,m+; 


1  /   n  n      v 

-     v  +■  v        )    . 

2  £,m+l  |,m'    ' 


(h.k.8] 


«/  i 


:+2»m 


5      1  5      1 

£+-ym+l  £+—  ,m-l 

~2Ah~ 
n  n  n  n 

l+Vm+l^^ _i.igLt-1  £+l,mrl  _"  J^m^l   , 

i4Ah 


(k.k.9) 


(5s>  1 

£,m+- 


£^m+l_     b£,m 
Ah 


(lj-.lMO) 


(Un 


y  .  i 


,n  ?n 

'£+l,m   _  £,m 
"Ah 


(lf.i^.ll) 


•5k- 


and 


,   vn  '   2 '      2     bl+l,m+l    l+l,m    £-l,m+l   "i-l,m 


(4.J+.12) 
are  already  available  from  previous  calculati 


Since  \k      and  ty      are  already  available  from  previous  calculations 


we  have 

,n     1  / ,n       ,n       ,n  ,n 

1   +  ^  1   +  ^    1  +  ^ 
--,m    i-g.m    2,n>f-    £,_ 


,n     1  / ,n       ,n       ,n       ,n 

i+^m    *-3>m    £>™4    ^m~^  )   •  (^-^.13) 


Similar  expressions  for  half  mesh  values  of  cp  can  also  be  derived.   As  for 
the  half  mesh  and  mesh  point  values  of  "Q  and  t,,   we  numerically  integrate 
Eqs.  (3.^.33)  -  (3.^.3^)  along  s  =  constant  and  get  two  simultaneous 
equations  in  t]      and  c. 


~n       ^h  /»  Nn   *-n       ^,n  ,,  /,  \n       Ah  ^n    /„  Nn 

i+p,m       J   '   ^+p^m    ^-p^m  £  >m        ^-pjm    i,m  ' 

and 

Ah  /.  \n    n       «,n  uX\                Ah  n    ,  „    >n 

■T  Vl,m\  1   +  K   1    =  ^  1   +"\  1   (?y}Mn  ' 

i+7j,m    U-,m  i--,m                i--,m     '       f4.lj-.15) 

With  the  values  of  T)   ,    and  £  given  by  previous  calculations,  the 

solution  of  Eqs.  (4.^.1^)  -  (k.k.Vi)  can  be  written 

,        g — 


2   y  «,">  (U.U.:6) 
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and 


..  /„  \n    n       L    r  Ah  ,„    >,n   -.21  «,n       Ah   ,,  Nn   ,„  Nn 


51 


i--,m    ^  '   ^   i«2,m 


'4-  i  +  [  f  <ty);  32 


(W.ifl 


i         /*  Nn      i+l,m    i-l,m  /.  ,  _,  n\ 

where     (£  ).   =  >  OAV *—  ,  (4.4.18) 

wy  i  ,m         2Ah  v      ' 


and 


(   }n     %m+l ^m-1  (4  4.19) 


We  also  have 


and 


n  n 

Slightly  more  complicated  expressions  for  rj      and  £      can  "be  derived 

I  ,m+2      *  ,nH"2 

by  using  the  same  method. 

The  set  of  equations  (4.4.1)  -  (4.4.21)  completely  describes  the 
numerical  iteration  procedure  for  two-dimensional  viscous  fluid  flow  problems. 
In  order  to  recover  the  transport  terms  of  each  flow  variable  for  asymptotical- 
ly steady  flow,  we  need  perform  the  "interpolation"  process.   In  our 
coordinates,  this  process  is  very  simple  as  expected.   Since  the  mesh  points 
are  all  moving  unidirectionally  from  left  to  right  along  curves  s  =  constant 
in  the  shock  reflection  problem,  it  is  only  necessary  to  use  formulas  like 
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(J+.2.7)  -  (^.2.1l)  to  interpolate  between  them.  The  boundary  conditions 
for  the  problem  shown  in  Figure  h.k  are  simply  as  follows,  thanks  to  our 
coordinate  system, 

(1)  Left  side  (uniform  incoming  flow) 

n      n        n      n        n      n       n      n 
P2,m  "  Pl,m   '   V2,m  =  Vl,m   '   ^2,m  '~~~   *l,m   '  P2,m  ==  Pl,m  ' 

VJ   =  0  ,  cp^   =0     for  m  =  1,2,..., M  ; 

(2)  Right  side  (continuous  outgoing  flow) 

n      n         n      n        „n     „n 

Pt     =  Pt   -,       )        Vt     =  Vt   -,     >         It     =  It  n     i 

L,m    L-l,m      L,m    L-l,m     L,m   -L-l,m 

n        n  *  n  o 

p    =  p  for  m  =  1,2,. . . ,M  ; 

L,m    L-l,m 

(3)  Top  side  (rigid  smooth  wall) 

n      n         n      n         n         n 
Pi,M  ==  Pi,M-l  '      Y£,M  =  Yi,M-l  '  *l,M   '~'~  *        *1,M-1  ' 


and 


n      n 


P*  m  =  P.  m_i     for  *  =  1,2,..., L  ; 


(*0   Bottom  side  (rigid  smooth  wall  with  straight  wedge) 

n      n       n      n       n         ,n 

Pi,l  =  PH,2  >      Vi,l  =  V*,2  '   *f,l  =  *  ~  *|,2  > 

p£,i  =  p?,2       for  i  =  1^2^--- >L  ; 


|--=Jt+2u-|     for  mesh  points  where  the  wedge  lies. 

:  appears  in  above  conditions  because  we  want  to  maintain  zero 
Lty  component  normal  to  the  walls.   In  actual  computation,  the  starting 
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half  wedge  angle  w  should  be  taken  small  enough  so  that  the  inclined  wedge 
face  does  not  intersect  any  of  the  initially  horizontal  streamlines.   It 
can  he  increased  gradually  after  every  few  time  cycles  to  achieve  the  desired 
final  angle. 

The  initial  flow  values  at  all  mesh  points  are  uniform  as  in 
Section  h.2.      For  simplicity,  we  take  them  to  he 

^°     in     °     -,  n    ,°     k  I    =   1,2,. ...L 

p*,m  =  1'°>   VMn  =  1-°>  h,^  =  2>      for  m  =  l,2;...;M  . 

The  value  for  initial  pressure  depends  upon  the  incident  shock  strength  | 
and  angle  a     specified  by  the  problem.   Using  oblique  shock  conditions  [10], 


we  have 


1  7  -  1*  z 1    ,  \ 

"i  = 7^~ '  :    (u-^> 


and 


?      vl.    7 

The  choice  of  v.    =1.0  implies  that  the  initial  speed  of  sound  CT  =  —  . 

X.    A  III  _!_      X  1— 

But  M  >  1  if  |   <  1  from  Eq.  (4.4.22).   Thus  we  are  guaranteed  to  have 
supersonic  incoming  flow.   The  desired  wedge  angle  w  for  generating  an 
incident  shock  with  angle  o^.  can  be  calculated  using  the  formula  [10], 

_L 

2    2 

M_   sin  a      -   1 

oa  =  tan"1  [  2  cot  a    g      ]  .  (4.4.21+) 

M  (7  +  cos  2a  )  +  2 

Other  constants  used  for  all  calculations  in  this  work  are 

7  -  1.4,  (j,  =  0.01,  L  =  50,  M  =  25,  Ah  =  0.05  and  At  =  0.02, 
where  At  has  been  chosen  to  satisfy  (4.3.2  ). 
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The  main  computer  program  is  written  in  FORTRAN  II  language  with 
inner  iteration  loop  coded  in  ILLIAC  II  machine  language.   We  test  for 
stationary  solution  of  a  given  shock  reflection  problem  by  comparing  all 
flow  variables  between  every  40  time  cycles.   If  the  relative  differences 
are  found  to  be  negligible,  we  stop  the  computation.   The  average  number  of 
time  cycles  required  for  stationary  solution  is  about  600  which  takes  about 
15  minutes  of  machine  time.   This  time  includes  starting  computation  from 
initially  uniform  flow  with  an  inclined  wedge  until  flow  becoming  stationary 
and  writing  flow  data  on  a  magnetic  tape  as  shown  in  Figure  4.5.   Figure  4.6 
gives  the  asymptotic  approach  of  stationary  solution.   The  changes  in  p  and 
v  along  streamlines  are  plotted  by  a  CalComp  67O  plotter. 

Main  results  of  this  work  are  summarized  in  Figures  4.7  -  4.9.   For 
given  incident  shock  strengths  £   =  0.9>  0.8  and  0.7,  we  plot  reflected 
shock  angle  cc-.   versus  given  incident  shock  angle  ax  obtained  by  our  numerical 
calculations.   We  compare  the  results  with  those  given  by  Von  Neumann's 
local  theory,  the  semi-global  theory  of  Section  2.3  and  Smith's  experiments. 
For  £   =  0.9,  we  also  give  numerical  results  for  the  case  of  no  explicit 
viscosity  effect,  i.e.,  u  =  0.   Numerical  results  that  we  obtain  by  the 
global  theory  are  in  much  better  agreement  with  Smith's  data  than  any  other 
results.   For  f   =  0.9,  the  results  given  by  the  1-Region  and  2-Region 
models  of  the  semi-global  theory  agree  quite  well  with  experimental  data. 
However,  the  agreement  becomes  worse  as  |   decreases,  i.e.,  as  the  incident 
shock  becomes  stronger. 

The  method  of  obtaining  reflected  shock  angles  from  our  computer 
is  demonstrated  in  Figure  4.5  for  the  case  |   =  0.9,  a     =   56  and 


>0.8  .   There  we  draw  constant  density  lines  on  a  density  variati 


on 
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chart  of  streamlines  (upper  diagram  of  Fig.  ^,5).   Transplanting  the  same 
constant  density  lines  onto  a  streamline  chart  in  physical  space  in 
rectangular  coordinates  (lower  diagram  of  Fig.  ^.5)>  they  become  curved. 
The  average  of  the  angles  of  constant  density  curves  made  with  the  wall 
measured  inside  the  reflected  shock  region  is  the  desired  reflected  shock 
angle.   The  same  technique  can  he  applied  to  velocity  and  pressure  variation 
charts  to  verify  the  results. 
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CHAPTER  V 
CONCLUSION 

The  purpose  of  this  work  was  to  show  that  the  Navier-Stokes 
equations  are  an  adequate  model  of  fluid  flow  in  the  case  of  the  weak  shock 
reflection  problem.   It  is  also  attempted  to  give  simpler  models  that  would 
agree  with  experimental  observations  and  be  computational  simple  so  that 
more  complex  problems  can  be  tackled. 

Previous  attempts  to  approximate  the  Navier-Stokes  model  in  various 
ways  have  failed  to  agree  with  experimental  results  for  weak  Mach  reflections, 
Many  of  the  efforts  were  discussed  in  Chapter  II.   One  of  the  first, 
Von  Neumann's,  which  uses  Rankine-Hugoniot  equations  locally  in  the  neighbor- 
hood of  the  shock  intersection  point  worked  for  strong  shock  reflections  and 
regular  reflections  but  not  for  weak  Mach  reflections.   This  is  presumably 
because,  in  the  presence  of  viscosity,  the  differential  equations  are  always 
elliptic  and  that  this  small  effect  becomes  more  noticeable  as  the  flow 
Mach  number  decreases  to  1.   Therefore  with  small  Mach  numbers,  the  effect 
of  the  downstream  boundaries  must  be  taken  into  account.   Many  authors  have 
made  modifications  of  this  local  theory  to  take  account  of  viscosity 
locally.   They  also  failed  to  work  for  weak  Mach  reflections,  presumably 
for  the  same  reason. 

In  Section  2.2„ ,  we  proposed  two  new  models  which  take  into  account 
the  shock  thickness  (due  to  viscosity)  and  the  effects  of  the  wall  to  which 
the  shock  is  attached.   Numerical  results  were  calculated  from  these  models 
for  a  range  of  incident  shock  angles  and  three  incident  shock  strengths. 
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While  these  were  pleasing  in  that  they  displayed  a  behaviour  similar  to 
Smith's  experimental  results,  they  were  not  as  close  to  those  results  as 
would  he  desired.   Furthermore,  the  more  complex  and  presumably  realistic 
model  gave  worse  results  in  this  respect  than  the  simpler  model.   However, 
this  so-called  semi-global  approach  appears  to  be  promising  in  that  an 
arbitrary  large  number  of  regions  can  be  used  to  represent  the  shock 
configurations,  but  at  this  time  a  better  understanding  of  the  flow  relations 
that  should  be  used  across  such  regions  is  needed. 

Since  these  models  were  based  on  simplifying  assumptions  made 
about  the  Navier-Stokes  model,  the  question  of  whether  or  not  this  is  an 
adequate  model  naturally  rises.   Previously  workers  in  this  field  have 
numerically  integrated  the  equations  of  motion  for  boundary  value  problems. 
They  either  ignored  viscosity  or  replaced  it  by  an  artificial  viscosity 
term  in  order  to  try  to  avoid  instability  that  often  results  near  shock 
regions.   Even  with  such  simplifying  assumptions,  the  computational  process 
is  quite  time-consuming.   Burstein  has  specifically  studied  the  shock 
reflection  problems.   His  approach  neglected  viscosity  (the  difference 
equations  serve  to  smear  out  the  shock  region  over  several  mesh  points 
much  as  viscosity  does  to  avoid  singularities).   His  results  agreed  well 
with  experiments  and  with  the  local  theory  for  strong  shock  reflections, 
but  not  for  weak  shock  reflections.   We  therefore  proceeded  on  the  assump- 

.on  that  the  true  viscosity  terms  should  be  retained.   Because  weak  shocks 
were  to  be  investigated,  it  seemed  reasonable  to  assume  that  the  viscosity 
;  was  constant  (the  temperature  does  not  change  much)  and  also 
heat  conduction  terms  can  be  neglected. 

[Terential  equations  were  numerically  integrated  in  a  mixed 
i.shioned  after  the  PIC  method  of  Harlow. 
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A  transformation  of  the  independent  variables  was  first  performed  in  order 
to  greatly  reduce  the  amount  of  computation  involved  in  this  two-stage 
process.   Numerical  results  were  obtained  for  a  range  of  incident  shock 
angles  over  a  set  of  three  incident  shock  strengths  whose  pressure  ratios 
were  0„9>  0-8  and  0„7  respectively.   This  computation  was  performed  on  the 
ILLIAC  II  computer  in  the  Department  of  Computer  Science,  University  of 
Illinois.   Each  shock  configuration  took  about  15  minutes  in  which  time 
about  600  time  steps  had  been  made  on  a  50  x  25  spatial  mesh.   The  inter- 
mediate and  final  flow  diagrams  were  plotted  by  the  computer  and  hand 
measurements  made  on  the  reflected  shock  angles.   They  agree  well  with 
Smith's  experimental  results  (see  Figs.  4.7  -  4.9).   Some  computations  were 
also  made  with  the  viscosity  terms  neglected.   The  results  did  not  agree 
well  with  experimental  data  (see  Fig.  4.7).   This  leads  us  to  suspect  that 
true  viscous  effect  is  indeed  important  for  weak  Mach  reflection  problems. 
It  also  shows  that  the  Navier-Stokes  description  of  fluid  flow  is  valid  in 
this  case. 

There  are  several  obvious  ways  to  improve  the  efficiency  of  the 
present  numerical  scheme.   For  shock  reflection  problems,  a  stationary 
incident  shock  could  be  maintained  in  the  two-dimensional  computational 
region.   This  is  because  the  generation  of  incident  shocks  by  placing  a 
wedge  in  a  uniform  flow  turned  out  to  be  quite  time-consuming.   As  for  the 
iteration  process  itself,  the  rate  of  convergence  may  be  increased  by  using 
appropriate  relaxation  schemes.   For  example,  in  treating  asymptotically 
stationary  problems,  the  time  is  no  longer  physically  meaningful  because  it 
is  used  in  iteration  to  avoid  numerical  difficulties.   We  can  consider 
variable  optimum  time  step  at  each  mesh  point  instead  of  a  uniform  time  step 
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for  the  whole  mesh.   This  local  time  step  may  serve  as  a  type  of  relaxation 
parameter  to  improve  the  convergence  rate.   Also  higher  order  difference 
approximations  may  be  used  to  reduce  the  number  of  mesh  points  required 
in  the  flow  region  and  hence  speed  up  the  computation. 
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